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LIGHT  SCATTERING  FROM  NONCONCENTRIC  SPHERES 


1  INTRODUCTION 

Atmospheric  aerosols  often  contain  small  inclusions  of  insoluble  material  which  can  affect 
their  scattering  properties.  These  inclusions  can  be  found  anywhere  within  the  host  aerosols. 
Chemical  and  biological  agents  can  also  exist  in  the  atmosphere  as  multi  component  aerosols. 
For  example,  bacterial  agents  can  be  dispersed  as  a  liquid  suspension  where  the  spore  or  cell  is 
embedded  in  a  water  droplet.  In  the  1950’s,  the  use  of  carrier  dusts  for  disseminating  chemical 
agents  was  tested.[l],[2]  Although  unsuccessful,  the  goal  of  this  work  was  to  modify  the  physical 
characteristics  of  the  chemical  agent  by  use  of  an  inert  dust  as  a  carrier  to  obtain  more  efficient 
aerosol  dissemination.  It  is  also  assumed  that  any  inclusion  or  encapsulation  can  directly  affect 
the  optical  properties  of  the  aerosol,  and  hence  its  scattering  characteristics.  This  influence  has 
consequences  for  remote  sensing  techniques  that  rely  on  the  aerosol  scattering  or  fluorescence 
for  detection  and  identification.  Therefore,  the  purpose  of  this  report  is  to  show  the  theoretical 
derivation  of  a  host  sphere  containing  one  spherical  inclusion  arbitrarily  located  within  the  host. 

Uniil  recently,  theoretical  predictions  of  scattering  for  a  host  sphere  containing  a  nonconccn- 

trically  positioned  smaller  sphere  was  not  feasible.  Most  experimentalists  had  to  be  content  with 
the  concentric  sphere  model, [3]— [5]  i.e.,  the  smaller  sphere  sharing  the  same  origin  as  the  host 
sphere.  Some  made  use  of  the  T-matrix  method  [6] -[8]  to  examine  a  radially  inhomogeneous 
sphere,  and  some  even  dabbled  with  the  effective  medium  approximations. [9]  It  was  not  until 
Friedman  and  Russek  [10]  published  the  addition  theorem  for  spherical  waves  that  the  solution 
to  the  nonconcentric  problem  was  even  realizable.  However,  Friedman  and  Russek  made  some 
errors  in  their  derivations.  Stein  [11]  (and  later  Cruzan  [12])  found  the  errors  and  came  up  with 
the  correct  translation  addition  theorem  for  spherical  wave  functions.  After  that,  the  solution  was 
attainable.  Fikioris  and  Uzunoglu  [13]  were  the  first  to  obtain  a  mathematical  expression  for  the 
scattering  of  a  nonconcentric  sphere  within  a  host  sphere,  followed  by  Borghese  et  a/.  [14]  Kirk 
Fuller  derived  the  relations  for  the  same  problem  using  an  order  of  scattering  approach.  [15]  In 
the  following  pages  we  derive  the  general  expressions  for  the  scattering  from  a  system  composed 
of  a  sphere  containing  a  nonconcentric  spherical  inclusion. 


2  THE  SOLUTION 

The  scattering  solution  of  a  sphere  containing  many  inclusions  is  difficult  to  solve,  and  even 
more  difficult  to  implement  because  of  the  limitation  of  present  day  computers.  For  this  reason, 
we  choose  a  simple  model  having  one  spherical  inclusion  arbitrarily  placed  inside  a  larger  sphere. 
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Figure  1 .  Geometry  of  the  nonconcentric  sphere  scattering  system.  A  host  sphere  of  radius 
fli  and  complex  refractive  index  mi  is  centered  on  the  xiyiZ\  coordinate  system.  The  inclusion 
sphere  of  radius  02  and  complex  refractive  index  m2  is  centered  on  the  X2y2Z2  coordinate  system 
at  a  position  xi  =  Q.  yi  —  0.  zi  —  d.  The  incident  radiation  is  a  plane  wave  traveling  in  the 
X  —  z  plane,  oriented  at  angle  a  with  respect  to  the  z  axis. 
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The  geometry  of  the  scatterer  is  shown  in  Figure  1 .  The  larger  sphere  (the  host)  of  radius  ci 
and  refractive  index  mi  completely  encloses  the  smaller  sphere  (the  inclusion)  of  radius  02  and 
refractive  index  m2-  The  center  of  the  host  sphere  is  situated  on  the  xi,  ?/i,  z\  coordinate  system. 
The  inclusion  sphere  is  centered  on  the  X2,y2-  22  coordinate  system  at  a  position  X]  =  0,  yi  =  0. 
zi  =  d.  Since  the  inclusion  sphere  is  centered  on  the  21  axis,  to  make  the  solution  completely 
general,  we  require  the  incident  radiation  to  impinge  upon  the  system  at  an  arbitrary  angle  a.  The 
problem  of  a  sphere  within  a  sphere  can  be  solved  by  simultaneously  satisfying  all  the  boundary 
conditions  at  the  surfaces  of  the  two  spheres. 

Figure  2  shows  the  physical  representation  of  the  fields.  Note  that  we  use  circular  arcs  to 
represent  spherical  traveling  waves  (spherical  Bessel  flmctions  of  the  3rd  and  4th  kinds)  and 
parallel  lines  to  represent  plane  waves.  The  boundary  conditions  at  the  surfaces  of  both  spheres 
are  satisfied  by  representing  the  fields  emanating  from  the  center  of  one  sphere  as  a  summation 
of  vector  harmonics  emanating  from  the  center  of  the  other  sphere.  We  use  the  vector  spherical 
harmonics  which  have  the  following  form: 


■  =  0 


im 
sin  B 


^3 


d  ^ 


N''^^  —  r- 


At, 


>)(Ar,)n(n  +  l)P„”'(cos^^)e*’"'^^ 


d  A, 


J__d 
kvj  drj 


'3 

im 


(1) 


where  the  index  j  corresponds  to  the  coordinate  system  used  (j  —  1,2)  and  z^^\krj)  are  the 
spherical  Bessel  functions  of  the  first,  second,  third,  or  fourth  kind  (p  =  1, 2. 3, 4).  and 


(2n+l)(n-m)! 
^  2{n  +  m)\ 


F-(cos^,)  = 

where  P^{cosBj)  are  the  associated  Legendre  polynomials. 


(2) 
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2.1  Host  Sphere 


We  now  examine  the  fields  which  strike  the  host  sphere  surface  (j  =  1).  We  consider 
an  arbitrary  field  incident  on  the  system  which  can  be  expanded  using  the  spherical  Bessel 
functions  of  the  first  kind,  j„(fcri).  We  will  examine  the  specific  cases  of  plane  waves  polarized 
perpendicular  and  parallel  to  the  y  axis  later.  The  incident  electric  field  may  be  expanded  as 

Einc  -  L  E  (3) 

n=0  m=— 77 

Similarly,  the  scattered  electric  field  may  be  expanded  using  the  spherical  Bessel  functions  of' 
the  third  kind.  hl^^{kri): 


Eic  =  i;  E  (4) 

n=0  777=— 77 

The  fields  near  the  boundary  \d\  <  vi  <  ai  inside  the  sphere  may  be  expanded  into  incoming 
and  outgoing  spherical  waves  using  spherical  Bessel  functions  of  the  fourth  kind  h!^\kiri)  and 
uiiru  kmu 


Ejph  “EE  J .  (5) 

n=0  77l=— 77 

The  application  of  boundary  conditions  at  the  host  sphere  surface  for  the  above  three  equations 
yields  tw'o  sets  of  equations; 


a^mhV''n{kai)  +  Cnmki^n\kai)  =  enmk^nHkidl)  +  «■!  ) '  (6) 

O.nm'H-'nikO-l)  +  CnmCn^H^“l)  =  CnmCn  )  +  5nmCn^H^l“l )  i  (”7) 

bnmWnikai)  +  dnm^n\kai)  =  H^l<^l):  (^) 

bnTnk\'^'niko.\)  +  dnmk\^n^\kai)  =  /nm^^n  \k\a\),  (9) 

where  -ipnir)  and  (9  =  1»2)  are  the  Riccati-Bessel  functions  defined  by 

fpr,  (r)  =  rj„  (r )  and  (r)  =  (r) ,  (1 0) 

and  the  primes  denote  derivatives  with  respect  to  the  argument. 
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2.2  Inclusion  Sphere 


Next  we  examine  the  fields  which  strike  the  inclusion  sphere  surface  (j  =  2).  The  fields 
inside  this  sphere  may  be  expressed  using  spherical  Bessel  functions  of  the  first  kind  jn{k7'r2)  '■ 

EL  =  f;  Y.  (11) 

n=0  fn=—n 

The  fields  near  the  boundary  outside  the  sphere  may  be  expressed  into  ineoming  and  outgoing 
spherical  waves  using  spherical  Bessel  functions  of  the  fourth  kind  h^^^{kir2)  and  third  kind 


Eext  —  ^  2  +  •SnmN^m.2  +  2  + 


(3) 


r(4) 


r(4) 


n=0  m=— 77 


Applying  boundary  conditions  at  the  sphere  surface  yields  two  sets  of  equations: 


(12) 


Pnmkl’ipni.k^Q'i)  —  f'nmk2^n  ^{^1(12)  +  ^{k\0,2)i 


I  Iiru  ^  n  \  - 


^riTTit^'^n  (^2^2)  ^nm^ri  ^(^1^2)  “f  (^1^2)' 

9nm^lt^n(^2*^2)  ~  ^nm^2‘sn  ''(^1*^2)- 


(13) 

(]4) 

(15) 

(16) 


We  can  eliminate  the  interior  field  coefficients  (pnm  and  qnm)  to  find  relationships  for  the  exterior 
field  coefficients.  After  a  little  bit  of  algebra,  we  have: 


_  ^  ki^'^'^>{kya2)ll!n{k2a2)  -  fc2^i‘'^H^TQ2)^/'n(^'2Q2) 

nm  k2^l^\kia2)tp'„{k2a2)  -  kl^n^\k■ia2)^pn{k2a2) 

_  ^  k2^',P^  (kia2)lpn{k2(l2)  ~  fcl^n^^(^TQ2)'0n(^'2Q2) 

ki^n  \kia2)ip'„{k2a2)  -  k2^n^\kia2)'ipn{k2a2) 


(17) 

-  ^77^717717 

(18) 

where  and  Q^,  are  the  Q  factors  which  contain  information  about  the  inclusion  sphere  such 
as  its  size  and  refractive  index. 


2.3  Fields  Interior  to  the  Host  Sphere 


The  fields  in  the  interior  of  the  host  sphere  are  expressed  by  equations  6-9,  while  the  fields 
exterior  to  the  inclusion  sphere  are  expressed  by  equations  17  £md  18.  At  first  glance,  it  is 
tempting  to  merely  equate  those  coefficients,  but  that  would  be  a  mistake  since  the  vector 
spherical  harmonics  expressed  by  those  equations  are  centered  about  two  different  coordinate 
systems.  This  is  where  the  importance  of  the  translational  addition  theorem  comes  into  play. 
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Stein[12]  and  Cnizan[13]  have  derived  translation  addition  theorems  for  vector  spherical  wave 
functions  which  can  be  used  to  express  the  coefficients  fnmi  Qnm^  STld  hfirrAT^  tCflllS  of  thc 
coefficients  Vj^rn^  ^nm'  inm^  u„rn-  Ajid  thc  Way  we  have  set  up  our  geometry  for  the  scattering 
system  now  comes  in  handy.  Because  we  require  that  the  inclusion  be  centered  along  the  z-axis, 
the  tw'o  origins  are  separated  only  by  a  displacement  in  z.  Hence,  there  is  no  rotation.  Thus,  for 
a  translation  along  the  z  axis  with  no  rotation,  the  vector  spherical  harmonics  are  related  by. 


M 

N 


(<?)  _ 
nm,2 

(19) 

Tj'—O 

(9)  _ 

n7Ti.2 

2^  ^n,n'  +  ^n.n' 

(20) 

77^=0 


where  q  denotes  the  order  of  the  spherical  Bessel  functions  {q  =  1,2. 3, 4).  This  relationship 
is  valid  in  the  region  where  r  >  |dl .  The  translation  coefficients  and  are  derived 

elsewhere.  [19]  We  note  that  the  expressions  are: 


.4 


(m.g) 


kid 


kid 
u  -r  i 


\ 


(77/-m+  l)(n/  +  m+  1)  ^^rn.g) 
y2ii'  -r  1  -r  oj 


(n'  -  m){ri'  +  m) 


n'  \\  {2n'  +  l)(2n'  -  1) 


^(Tn,gj 


(21) 


(rn,g)  _  -jki'md  f^jm.g)  (22) 

-^77  77^  //  /  i  “iN  Tt.T)^ 

n'(77^  +  lj 

The  are  scalar  translation  coefficients.  Recursion  expressions  for  these  scalar  coefficients 

can  be"derived  using  the  method  of  Bobbert  and  Vlieger.[17]  The  details  are  shown  elsewhere.[19] 
The  necessary  scalar  translation  coefficients  are: 


-  ^/2^^^3Akldl 
=  -^2n'  +  lj„~{hd)- 


M0,g) 


2n  4"  3 


(^n)V2n'  +  l 


2n  +  1 


y(0,g:) 


+ 


n\ 


I2n'  +  1 
2n 


l/o(o.g)  I  1^  /2n  +  1  ^(o.g)  1 

-  ("  +  lly  2n'  +  3^”-"'+^  J  ■ 


(23) 

(24) 


(25) 
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y'(r?  -  777  +  !)(??  +  m)(277'  +  =  ■y/(r?'  -  m  +  l)(n'  +  m){2n'  +  - 

j  in'  -  rn  +  2)(r)'  -  rn  + 

'  \  (2n'  +  3) 


(n'  +  m)(r7'  +  777  -  1) 

(STTT) - ■  P6) 


'^71. 7?^  '^72.7?' 


From  these  equations,  we  see  that 


_  Am, 4}  _  >(-m,3)  _  .(m) 

"'"^71. rz'  ^72, n'  ^7z,n'  ’ 


p(7r7.3j  _  j^imA)  _  p(-77i,3j  _  p(77i) 

^n.n'  ^n.n'  ^n.n' 

^(771,3)  _  ^{ttiA)  _  ^(-m,3)  _  ^{m) 

^72.7?'  ~  ^n,7?'  “  ^71,72'  “  ^77.72'* 


The  advantage  of  expanding  the  interior  fields  of  the  host  sphere  in  the  third  and  fourth 
kinds  of  spherical  Bessel  functions  rather  than  the  first  and  second  kinds  is  that  only  one  set 
of  translation  coefficients  need  to  be  calculated.  Using  equations  19,  20,  and  substituting  those 
expressions  into  equation  12,  we  obtain 


=  E  E 


n=0  m  =  ~n 


(772,3) -VT(3.j 


'^^mn  2-^  ^71.72'  ^72,72' 

_72'  =  0 

-1-7  ^  -A- 

+tnm  2_^  1  +  £>„_„/ 


+  u„„  2:  +  4"- 'n2; 


The  field  outside  of  the  inclusion  sphere  is  now  expressed  in  terms  of  the  vector  spherical 
harmonics  of  the  host  sphere.  Therefore,  it  is  now  possible  to  equate  equation  29  with  equation 
5.  Looking  only  at  the  we  see  that 


E  E  ^nrM^  =  E  E  E  Kn'  +  ^nm  E 


72=0  m——n 


r2=0  772  =  -r2  L  72'=0 
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Multiplying  both  sides  of  this  equation  by  f  and  using  the  orthogonality  condition 

{^ni^mj  on  the  left  hand  side,  and  on  the  right  hand  side),  we  obtain 


DC 

“I" 

77  =  0 


which  can  be  rewritten  as 

oc 

E  Am  i  Q 

*Ti'm^n' ,71 

77 '=0 

Continuing  the  same  line  of  reasoning,  we  can  show 


fn 


77'=0 


4777 

n'm^n^n 


(31) 


(32) 


(33) 


DC 

9nm  ~  'y  ]  ^Tt'm-^n'.n  "f  ’ 

n'=0 


(34) 


and 


oc 


.b:- 


n5> 


n'=0 

These  are  the  four  magic  equations  which  relate  the  coefficients  in  the  two  different  coordinate 
systems.  Armed  with  these  new  expressions,  it  is  now  |x>ssible  to  replace  one  set  of  coefficients 


with  another. 


Substituting  equations  32  and  34  into  equation  6  we  have 

OC 

ClnmklWn{kcil)  +  ^  + 


Tl'=0 

cx: 


tn'mA^Zn'  +  (36) 

T7'=0 


Using  equations  17  and  18  for  and  Sn'm.  we  can  substitute  those  expressions  into  the  above 
equation  to  obtain 


^nm  V;„(A:ai)  +  CnmCn 


n'=0 

[«f  (fca.)  +  ■  (37) 


Likewise,  we  can  substitute  equations  32  -  35  into  equations  7-9  (and  together  with  equations 
17  and  18)  to  obtain  the  following  set  of  simultaneous  equations,  which  can  be  solved  to  yield 
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the  scattering  coefficients  (0^,^  and  d^rn)  or  the  internal  field  coefficients  and  ?/„„,)  of  the 
host  sphere: 


^nm  ^‘77  (A'ai)  +  Cnm^'^^^{kai) 


^nm  V-^n  (koi)  -f  d  nm  ^n’ikai) 


^nm  Wji  (koi)  +  d 

nm 


7?'=0 


n'=0 

tl'Tl 


+ 


*'1  „'=0 

7/  /  4^ 

^71  m  ^71,1 


(38) 


(39) 


(40) 


2.4  Scattering  Coefficients 


The  set  of  four  equations  (37  -  40)  is  the  final  product  in  our  search  for  the  scattering 
coefficients.  We  have  four  equations  and  four  unknowns,  the  unknowns  being  Cnrn^  dnm-  inm- 
and  Unm-  The  interior  field  coefficients  inside  the  host  sphere,  tnm  and  Unm,  (which  can  also  be 
viewed  as  the  exterior  field  coefficients  of  the  inner  sphere)  may  be  calculated  by  eliminating 
the  scattering  coefficients  c„m  and  dnm  in  equations  37  -  40.  Since  all  the  elements  comprising 
the  matrix  in  the  solution  are  known  except  for  the  coefficients  and  Unm-  we  can  perform 
an  LU-decomposition  to  solve  for  those  interior  field  coefficients.  The  matrix  representation  is: 


where 


_  .  rp{Tn,l}  ,  Tr(m,lJ 

^nmln  /  ^  n,n'  '  ^n'm^n,n'  ’ 

n'=0 

(41) 

DC 

h  —  i  rp{m,2)  jj{rti,2} 

n^=0 

(42) 

(43) 

AZn’  {kii'Hka,)  [?f(i-.a,) 

(44) 
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(45) 


BZn'  {kiC^'Hkar)  [eHhai)  + 


TT(m.l) 

^n.n' 


Bn,n'  [^fH^lOl)  +  Qn'^^n\^l(^l)] 

-he^ika,)  [ef '(fciai)  +  Hfciai)]  }  • 


K.n'  {h^n^Hkai) 

[e\k,a,)  +  • 


(46) 


(47) 


The  scattering  coefficients,  c„m  and  d„m  may  then  be  calculated  using  equations  37  -  40. 
Once  we  know  the  scattering  coefficients,  we  can  determine  the  scattering  matrix  elements  and 
other  parameters  associated  with  the  scatter. 

2.5  Plane  Wave  Expansion 


Since  tl'ic  small  sphere  is  centered  on  the  z  axis,  for  the  solution  to  be  completely  general, 

the  incident  plane  wave  must  be  in  an  arbitrary  direction  in  the  x-z  plane.  Two  plane  waves  must 
be  considered  in  order  to  account  for  each  polarization  state.  When  the  plane  wave  is  polarized 
perpendicular  to  the  x-z  plane  (TE),  the  coefficients  are  found  to  be  [18]; 


= 

^nm 


^J{n  -  m){n  -t-  m  +  l)P^^\cosa) 


n(r}  +  1)  L 
yj [n  -  m  +  1)(t?  +  m) ~ ^ (cos a ) 

27”+2  Q  _ 

-P;^{cosa). 


n{n  +  1)  da 


(48) 

(49) 


b 


nm 


^TE 

nm 


f"+^(2n  +  1) 
n(n  -f  1) 


-m  +  l){n  -m  +  2) 
(2n  +  l)(2n-h3) 


(cosq) 


(n  -H  m  -h  l)(n  -H  m  -f-  2) 
(2n  +  l  )(2n  -h  3) 


2^^“^  mP”‘(cosQ) 
n(T? -f  1)  sino 


(50) 

(51) 
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When  the  plane  wave  is  polarized  in  the  x-z  plane  (TM),  the  coefficients  are  found  to  be 


^nm  ^nTT7 


bnm  = 


2.6  The  Scattering  Amplitudes  and  Efficiencies 

We  consider  the  scattering  amplitudes  in  the  far  field,  where  kri  »  ka.  The  scattered  fields 
in  this  case  are  in  the  6  and  (p  directions.  In  this  limit,  the  spherical  Hankel  function  reduces  to 
spherical  waves: 

h^^^{kT)  ~  (54) 


The  scattering  amplitudes  can  be  expressed  in  the  form  of  the  matrix, 
/  eT  \  /  5i  54  \  (  ETe  \ 


—ikri  \  Sz  S2  )  \  Ej 


The  scattering  amplimde  matrix  elements  are  solved  by  expanding  the  scattered  electric  fields 
(equation  4)  in  terms  of  the  vector  wave  functions  and  then  expanding  the  vector  wave  functions 
(equation  1 )  in  terms  of  the  polarization  directions.  Since  we  are  concerned  with  the  far  field, 
we  may  drop  the  f  component  and  keep  only  the  9  and  (p  components.  After  some  algebra,  we 

oc  «  r  Q  ' 

s,  =  f:  f;  X  +  .  (56) 


5.  =  f:  d^^£^P-(cos9,}  +  c^^^P-{cos9,)  . 

«=-'T;  t  X  c™~P„"(cos9,)+C,'^C(™s«.) 

=  i:  (-ire'”'-' X  + 

n=0m=“r7  L  ^  ^ 


f  ±  X  UM^P»(cosft)  +  c™| 

n=0  Tn=—n  L  ^ 


7^=U  Tn  =  —Ti  L  ^  J 

Using  the  following  relationship  for  normalized,  associated  Legendre  polynomials, 

F-”*(cos0i)  =  (-l)’"4"^(cos0i), 

the  following  relationships  between  the  scattering  coefficients  may  be  derived: 

PTE  _  /  ■,\mJrE 


^nm  \  ■*'/  ^nTTi  ■ 
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=  ( 

_-j  \7774-l^rA/ 
■*■  /  ^  nm 

77777,  ' 

=  ( 

11 

(-irC- 

(61) 


where  m  =  -m. 


The  scanering.  extinction,  and  absorption  efficiencies  of  the  system  are  defined  as  the  cross 
sections  per  projected  area  and  may  be  expressed  as 


Q  sea 


Qcxi 


_9 


{ka,y 


(A’a,)2 


X  Re 


X 


.77  =  1 


+  i;  (|c;:.sr+  c™|-  + 


jTM 

^nm 


(62) 


/  ,  1  \  TE^TE*  i  jTElTE*  ,  ^TM ^TM*  , 

T\\tI  +  1)  /  ^  {^nm^nm  ^nm^nm  ^nm  ^nm  ^nm  nm  j 


ln=l 


lit IL  ^  ^  ^ 

(63) 

Qahs^Qexi-Qsca-  (64) 

where  a;„  and  are  the  complex  conjugates  of  a„„  and  Km,  respectively.  The  asymmeby' 
parameter,  g,  can  be  expressed  as 

4 


9  = 


Qscaika)'^  m 
+n(r7  +  2), 


(J'EjTE*  ,  ^TMjTMA 

tnRe  4-  d^^  j 


(n  —  m  +  l)(r?  +  m  +  1) 


X 


(2n  +  3)(2n  +  l) 

\xl„TE^TE*  J_  jTExTE^  ,  TM  TM*  ,  jTMjTM*  1 


(65) 

Detailed  derivations  for  the  asymmetry  parameter  and  the  efficiencies  are  given  elsewhere. [19] 


3  SPECIAL  CASES 


Under  certain  conditions,  the  equations  describing  the  internal  and  scattered  fields  may  be 
simplified.  In  this  section  we  examine  some  of  these  cases. 

3.1  Simplifications  of  Ql  and 

When  the  inclusion  satisfies  eertain  conditions,  equations  17  and  18,  which  describe  the  co¬ 
efficients  and  Q® ,  can  be  simplified.  In  this  subsection  we  demonstrate  these  simplifications. 
First,  when  the  optical  size  of  the  inclusion  is  much  smaller  than  the  wavelength  (^202  <s:  1). 

(66) 
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^  ,  4?(A:2Q2)^  (nil-rm) 

”  3  (ml  +  2mi) 


(67) 


These  are  similar  to  the  scattering  coefficients  for  a  Rayleigh  scatterer.  If  the  inclusion  is  a 
perfectly  conducting  medium  (m2  — >  00).  the  coefficients  can  be  reduced  to 


II 

^\kia2) 

(68) 

Qn  = 

ef^(Ayc2)' 

(69) 

When  the  perfectly  conducting  inclusion  is  much  smaller  than  the  wavelength  (ao  <  A), 
equations  68  and  69  further  reduce  to 


2i(kia2Y 

3  ■ 

(70) 

^  4i(A:ia2)^ 

(71) 

3.2  Mie  Reduction 


It  should  be  noted  that  if  the  refractive  index  of  the  inner  sphere  is  equal  to  that  of  the  outer 
sphere,  then  the  exterior  field  coefficients  of  the  inner  sphere,  and  are  umty  (equations 

1 7  and  1 8).  The  Riccati-Hankel  functions  of  the  first  and  second  kind  on  the  right  hand  sides 
of  equations  37  -  40  combine  to  form  Riccati-Bessel  functions  of  the  first  kind: 


.  2k  ^ 

Onm'^nikO])  j4„  „/'i/^n(A'l  til  )  -|-  ) . 

^’1  n'=0 


(72) 


On7n^'n(A:Oi  )  -|-  (kOl)  —  2  ^  -|-  Uji'mBn.n'Vnikiai)  ■ 


n'=0 

DC 


bnm'kyn{ko.\)  -f  —  2  ^  Un'Tn-^n,n''^n{k\ai)^ 


r)'=0 


2k  ^ 

brun'fp'nika-l)  +  dnm^n^Kkai)  =  —  ^  tn'mB'^n''^nik\(l\)  +  ^^n'm^n.n'V'n(A^l°l )  • 

n'=0 


(73) 

(74) 

(75) 


Physically,  the  Riccati-Hankel  functions  of  the  first  and  second  kind  represent  outgoing  and 
incoming  spherical  waves,  respectively.  With  the  inner  sphere  removed,  there  is  no  scattering 
object  within  the  larger  sphere,  and  standing  waves  represented  by  Riccati-Bessel  functions  of 
the  first  kind  are  contained  within  its  interior.  Note  that  substituting  equations  34  and  35  into 
the  right  hand  side  of  equations  72  -  75  yields  the  boundary  equations  for  Mie  theory. 
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3.3  Concentric  Spheres 


WTien  the  host  and  inclusion  spheres  are  concentric,  i.e.  they  share  the  same  origin  {d  =  0), 
the  translation  coefficients  are  simplified,  and 

BZ-  =  (76) 

where  bnn'  is  the  Kronecker  delta  function.  The  scattering  coefficients  described  by  equations  37 
-  40  now  reduce  to  the  accepted  concentric  spheres  expressions  given  by  Aden  and  Kerker.[4] 


4  RESULTS  AND  DISCUSSIONS 

In  this  section  we  examine  the  scattering  of  the  nonconcentric  sphere  system  as  a  function  of 
system  parameters.  As  in  any  new  theory,  confidence  of  a  new  result  relies  on  its  simplification  to 
previously  known  solutions.  But  rather  than  analytically  simplifying  the  general  nonconcentric 
expressions  into  specialized  cases  (we  have  already  shown  the  special  cases  in  the  previous 
section),  we  will  keep  eveiy^thing  general  and  allow  the  computer  to  numerically  give  us  the 
desired  results.  The  greatest  advantage  of  doing  it  this  way  is  the  assurance  that  the  program 
is  coded  correctly;  or  at  the  very  least,  the  code  reduces  correctly  in  the  specialized  cases. 
Simplifications  of  the  system  can  be  made  by  assigning  to  it  special  values  so  that  comparisons 
with  Mie  theory  and  concentric  spheres  theory  may  be  done.  We  will  show  that  when  we  simplify 
the  system,  we  obtain  numerical  values  that  are  identical  to  Mie  theory  and  concentric  spheres 
theory.  We  then  consider  two  specific  but  arbitrary  cases  of  the  noncencentric  system.  The  first 
is  an  air  inclusion  (m2  =  1.0  +  Oi)  in  a  host  sphere  composed  of  water  (rbi  —  1.33  +  Of).  The 
second  is  a  perfectly  conducting  inclusion  in  a  host  sphere  composed  of  water.  In  both  cases  the 
host  radius  is  equal  to  the  wavelength  of  the  incoming  radiation  (a-i  =  A). 

4.1  Comparison  with  Mie  Theory 

When  the  refractive  index  of  the  host  (mi)  and  the  refractive  index  of  the  inclusion  (m2) 
are  identical,  we  can  say  that  the  system  is  homogeneous.  In  this  case,  the  nonconcentric  sphere 
should  reduce  to  Mie  theory.  We  take  the  host  and  the  inclusion  to  be  glycerol  spheres  with 
refractive  index  m  =  1.4746  +  Oi  at  wavelength  A  =  0.5145//.m.  Figure  3a  shows  the  scattering 
intensity  of  the  glycerol  sphere  as  a  function  of  the  sphere’s  radius.  This  plot  is  obtained  using 
the  Mie  code  presented  by  Bohren  and  Huffman.  [20]  The  familiar  quasi-periodic  resonance 
peaks  are  readily  observable.  Figure  3b  shows  the  same  plot,  but  using  the  nonconcentric  sphere 
model.  We  note  that  all  the  cycles  of  resonances  are  the  same;  in  fact,  we  have  normalized  the 
scattering  intensity  so  that  numeric  values  from  the  Mie  code  and  the  nonconcentric  code  are 
the  same  (as  can  be  seen  in  Figure  3).  Also  listed  in  Table  1  are  the  results  from  a  run  using  the 
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Radius()j,m) 


Figure  3.  Comparison  of  the  scattering  intensity  as  a  function  of  droplet  radius  (m 
1.4746  +  Or)  using  the  Mie  code  from  Bohren  and  Hufftnan  (a),  and  the  nonconcentric  code  (b) 
with  7/1 1  =  mo  =  1.4746  +  Of  The  resonance  structures  are  cleariy  evident  in  both  cases. 


TABLE  1 


Comparison  with  Mie  Scattering  output 


Input: 

REFRE  =  1 .550  REFIM  =  0.0 

SPHERE  RADIUS  =  0.525  pm 
WAVELENGTH  =  0.6328  pm 


Output: 


Osca  =  3.10543  Qext  =  3.10543  g  =  0.633137 


angle 

Sll 

POL  (SI 2) 

S33 

S34 

0.0 

1.00000 

0.0 

1.00000 

0.0 

9.0 

0.785391 

4.5981  lE-03 

0.999400 

3.43261  E-02 

18.0 

0.356897 

4.58540E-02 

0.986022 

0.160184 

27.0 

7.661 19E-02 

0.364744 

0.843602 

0.394077 

36.0 

3.55355E-02 

0.534997 

0.686967 

-0.491787 

45.0 

7.01845E-02 

-9.59953E-03 

0.959825 

-0.280434 

54.0 

5.74324E-02 

-4.77927E-02 

0.985371 

0.163584 

63. U 

2.19660E-02 

0.440604 

0.648043 

0.621216 

72.0 

1.25959E-02 

0.831996 

0.203255 

-0.516208 

81.0 

1.73750E-02 

-3.41670E-02 

0.795354 

-0.605182 

90.0 

1.24601E-02 

-0.230462 

0.937497 

0.260742 

99.0 

6.79093E-03 

0.713472 

-7.17406E-03 

0.700647 

108.0 

9.54240E-03 

0.756255 

-3.94746E-02 

-0.653085 

117.0 

8.634 19E-03 

0.281215 

0.536251 

-0.795835 

126.0 

2.2742  lE-03 

0.239612 

0.967602 

7.95797E-02 

135.0 

5.43998E-03 

0.850803 

0.187531 

-0.490882 

144.0 

1.60244E-02 

0.706334 

0.495255 

-0.505781 

153.0 

1.88853E-02 

0.891081 

0.453277 

-2.268 17E-02 

162.0 

1.95254E-02 

0.783320 

-0.391613 

0.482752 

171.0 

3.01676E-02 

0.196194 

-0.962069 

0.189556 

180.0 

3.83189E-02 

0.0 

-1.00000 

0.0 
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TABLE  2 

Comparison  of  Nonconcentric 
with  Coated  Sphere  Theory 


Input: 

X  =  3.0pm,  mi  =  1.409  +  0.1747i,  m2  =  1.59  +  0.66i 


Output: 


Radius  (urn) 

Coated  Sphere 

Nonconcentric 

a,  =  6.2650 

Qext  =  2.32803 

Qext  =  2.32803 

a,  =  0.1710 

Qsca=  1.14341 

(ica=  1.143413 

a,  =  5.50 

Qex.  =  2.34481 

Qext  =  2.344814 

a2  =  1.50 

Qsca=  1.13227 

Q5ca=  1.132269 

ai  =  3.50 

Qext  =  2.58325 

Qex,  =  2.58325 

a.  =  2.00 

Qsca  =  1.28485 

(iea  =  1.28485 

Input: 

X  =  3.0pm,  mi 

=  1.409  +  O.Oi,  m2  =  1.592  +  0.66i 

Output: 

Radius  (pmi 

Coatea  Sphere 

Nonconcentric 

a,  =  6.2650 

Qext  =  2.77588 

Qext  =  2.77588 

a,  =  0.1710 

Qsca  =  2.775 18 

Qsca=  1.77518 

a,  =  5.50 

Qext  =  2.04387 

Qext  =  2.04387 

a,  =  1.50 

Qsca  =  1.85917 

Q,ca  =  1.859173 

a  =  3.50 

Qex,  =  3.19478 

Q,,,  =  3.194786 

a.  =  2.00 

Q,ca  =  2.41 198 

Q,ea  =  2.41 1986 

Input: 

X  =  3.0pm,  m. 

=  1.409  +  O.Oi,  m2  =  1.592  +  O.Oi 

Output: 

Radius  (um) 

Coated  Sphere 

Nonconcentric 

a,  =  6.2650 

Qext  =  2.77604 

Qext  =  2.77604 

a,  =  0.1710 

Qsca  =  2.77604 

aea  =  2.77604 

ai  =  5.50 

Qext  =  2.31972 

Qext  =  2.319716 

a,  =  1.50 

(ica  =  2.31972 

Qica  =  2.319716 

a,  =  3.50 

Qext  =  2.28173 

Qext  =  2.28173 

a:  =  2.00 

Qsca  =  2.28173 

aca  =  2.28173 
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nonconcentric  code  with  inputs  taken  from  page  482  of  Bohren  and  Huffman.[20]  The  numerical 
values  for  both  programs  are  identical.  This  tells  us  that  our  nonconcentric  sphere  model  does 
indeed  reduce  to  Mie  theory  in  the  appropriate  limit. 

4.2  Comparison  with  Concentric  Theory 

When  the  displacement  d  separating  the  two  origins  is  zero,  we  have  a  concentric  sphere 
system.  Using  the  concentric  code  presented  in  Bohren  and  Huffman,  we  are  able  to  show  that, 
in  all  cases,  our  nonconcentric  sphere  model  reduces  to  the  concentric  sphere  model.  Sample 
results  of  using  the  two  methods  are  given  in  Table  2  (the  first  set  of  output  can  be  compared 
to  page  489  of  Bohren  and  Huffman). [20] 

4.3  Examination  of  Two  Types  of  Inclusion 

In  this  section,  we  examine  the  scattering  properties  of  the  nonconeentrie  sphere  model  by- 
considering  two  specific  but  arbitrary  cases.  The  first  is  an  air  inclusion  (m2  =  1.0  +  O.Oi)  in 
a  host  sphere  composed  of  water  (mj  =  1.33  +  O.Oi).  The  second  is  a  perfectly  conducting 
inclusion  in  a  host  sphere  composed  of  water.  In  both  of  these  cases,  the  host  radius  is  equal  to 
the  wavelength  of  the  incident  radiation  (ci  =  A). 

Figure  4  shows  the  extinction  efficiencies  as  a  function  of  displacement  distance  d  when 
the  incident  angle  a  =  0°  and  the  beam  travels  parallel  to  the  z  axis  for  inclusion  spheres  of 
different  radii.  For  comparison,  we  also  plot  the  extinction  efficiencies  for  a  sphere  composed 
entirely  of  water  (0,2  =  0,  shown  by  the  solid  lines,  top  and  bottom),  for  which  Qext  —  3.916. 
For  a  sphere  composed  entirely  of  a  perfectly  conducting  medium  (top  graph),  Qext  =  2.094.  For 
a  sphere  composed  entirely  of  air,  we  have,  of  course,  Qext  —  0.0.  The  most  remarkable  feature 
of  the  graphs  is  their  symmetry.  The  extinction  is  independent  of  the  sign  of  the  displacement; 
i.e.  Qcxt[(^}  ~  Qexti 

The  result  that  Qext{d)  =  Qext{-d)  is  a  direct  consequence  of  the  reciprocity  theorem  which 
states  that  if  a  scattering  system  whose  scattering  matrix  is  given  by  equation  55  is  replaced 
by  its  mirror  image  (the  mirror  is  perpendicular  to  the  direction  of  the  incident  field),  then  the 
scattering  due  to  the  mirror  system  in  the  forward  direction  {61  =  a,  (^1  =  0°)  is  given  by 

F;-(a,o^)  \  _  ^  (  s,  -53  W  \  .77. 

£:r(«,0O)  j  1,  -54  S2  )\ErM)' 

We  know  from  the  optical  theorem  that  the  extinction  is  proportional  to  the  real  part  of 
the  scattering  amplitude  in  the  forward  direction,  and  from  equations  48-53  and  56-61,  it  can 
be  shown  that  there  is  no  coupling  of  modes  in  the  scattering  plane,  i.e.,  S3  and  do  not 
contribute  to  the  scattering  amplitude  when  (^  =  0*^.  Therefore,  the  extinction  of  the  scattering 
system  discussed  in  this  paper  is  equal  to  the  extinction  of  its  mirror  image. 
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d(X) 

Figure  4.  The  extinction  efficiencies  as  a  function  of  displacement  distance  d  when  the 
incident  angle  q  =  0°  for  a  perfectly  conducting  inclusion  (top)  and  an  air  inclusion  (bottom) 
inside  an  ai  =  A  water  host  (mi  =  1.33  +  Oz)- 
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Figure  5.  The  backscattering  total  intensity  5ii  as  a  function  of  displacement  distance  d 
when  the  incident  angle  a  —  0""  for  a  perfectly  conducting  inclusion  (top)  and  an  air  inclusion 
(bottom)  inside  an  O]  =  A  water  host  (mj  =  1.33  -h  Oz). 


Figure  5  shows  the  backscattering  intensities  (incident  angle  0  =  0*^  and  6i  —  180*^)  as  a 
function  of  displacement  distance  d  for  the  same  cases  shown  in  Figure  4.  Also  shown  on  both 
graphs  are  the  backscattering  intensities  when  there  is  no  inclusion  (02  =  0),  in  which  case 
the  backscattering  5ii  =  1.837.  Unlike  the  plots  shown  in  Figure  4,  these  plots  do  not  show 
a  symmetric  Su.  There  are  also  some  differences  in  the  backscattered  intensities  between  the 
two  systems  in  this  figure.  For  the  case  of  the  perfectly  conducting  inclusion,  the  backscatter 
is  generally  greater  than  for  a  sphere  with  no  inclusion  (at  times  by  more  than  an  order  of 
magnitude);  whereas,  for  the  case  of  the  air  inclusion,  the  backscatter  generally  stays  within  a 
factor  of  two  of  the  backscatter  of  the  system  with  no  inclusion. 

The  difference  in  the  backscattering  between  a  system  with  a  perfectly  conducting  inclusion 
and  an  air  inclusion  can  be  bener  understood  by  considering  the  rays  which  strike  the  system. 
Incident  rays  striking  the  surface  of  the  host  sphere  either  refract  convei^ently  within  the  sphere 
(for  our  case  when  the  refractive  index  of  the  host  sphere  is  greater  than  the  refractive  index  of 
the  incident  medium)  or  reflect  off  the  outer  surface.  Some  of  the  refracted  rays  which  reflect  off 
the  surface  of  the  inclusion  may  be  traced  into  the  backscatter  direction  (these  are  the  rays  which 
strike  the  inclusion  at  normal  incident  or  strike  the  inclusion  at  2:2  =  ya  =  0).  The  reflectance 
from  the  perfectly  conducting  inclusion  is  100%;  whereas,  the  reflectance  from  an  air-water 
interface  is  approximately  2%  near  normal  incidence.  Because  of  the  greater  reflectance  at  the 
conductor  interlace,  we  expect  the  backscatter  to  be  greater  when  there  is  a  perlectly  conducting 
inclusion.  The  reflectance  at  the  air  inclusion  interface  is  approximately  the  same  as  at  the  outer 
surface  interface,  so  we  do  not  expect  the  backscatter  for  the  air-inclusion  system  to  var\'  as 
greatly  from  the  sphere  system  without  the  inclusion. 

Figure  6  shows  the  extinction  efficiencies  for  an  02  =  A/4  inclusion  as  a  function  of  the 
incident  angle  a  for  four  displacements.  When  the  inclusion  and  host  spheres  are  concentric,  the 
extinction  does  not  vaiy',  as  expected  from  symmetry.  As  the  inclusion  distance  increases,  the 
magnitude  of  the  total  variation  in  the  extinction  becomes  greater.  Also,  the  number  of  inflection 
points  in  the  extinction  curves  increases.  Essentially,  the  closer  the  inclusion  is  to  the  host  sphere 
surface,  the  greater  the  effect  it  has  on  the  extinction  of  the  system.  This  point  is  reinforced  in 
the  next  section  where  we  examine  the  resonance  of  the  composite  system.  We  also  note  that 
because  of  the  system  symmetry,  the  extinction  efficiencies  must  be  symmetric  about  o  =  0°. 
and  from  the  reciprocity  theorem,  they  must  also  be  symmetric  about  q  =  90°. 

Figure  7  shows  the  backscattering  intensities  (0i  =  tt  -F  q)  for  an  02  =  A/4  inclusion  as 
a  function  of  incident  angle  q  for  the  same  displacements  as  in  the  previous  figure.  As  was 
the  case  for  the  extinction  efficiencies,  the  number  of  inflection  points  in  these  backscattering 
intensity  curves  increases  as  the  inclusion  distance  increases.  However,  the  symmetry  about 
a  =  90°  which  exists  in  the  extinction  efficiencies  of  Figure  6  is  not  present  in  the  backscattering 
intensities. 

Figure  8  shows  the  extinction  efficiencies  as  a  function  of  the  inclusion  radius  for  four 
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Figure  6.  The  extinction  efficiencies  as  a  function  of  the  incident  angle  q  for  a  perfectly 
conducting  inclusion  (top)  and  an  air  inclusion  (bottom)  inside  an  O]  =  A  water  host  (rh]  = 
1.33  +  Oz).  The  inclusion  radius  is  a2  =  A/4. 
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Figure  7.  The  backscattering  total  intensity  Su  as  a  function  of  Ae  incident 
perfectly  conducting  inclusion  (top)  and  an  air  inclusion  (bottom)  inside  an  Ci  — 
(777^  =  1.33  4-  0?).  The  inclusion  radius  is  oo  =  A/4. 
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Figure  8.  The  extinction  efficiencies  as  a  function  of  inclusion  radius  when  the  incident 
angle  q  =  0°  for  a  perfectly  conducting  inclusion  (top)  and  an  air  inclusion  (bottom)  inside  an 
Q]  =  A  water  host  (rh)  =  1.33  +  Or). 
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Figure  9.  The  extinction  efficiencies  as  a  function  of  inclusion  radius  when  the  incident 
angle  q  =  0°  and  aj  =  q2-\-  d  for  a  perfectly  conducting  inclusion  (top)  and  an  air  inclusion 
(bottom)  inside  an  Oi  =  A  water  host  (mj  =  1.33  +  Oz)* 


different  separation  distances  when  0  =  0*^.  The  solid  lines  in  both  graphs  (d  =  0)  correspond  to 
the  cases  where  the  host  and  inclusion  spheres  are  concentric.  We  also  calculated  the  extinction 
efficiencies  for  systems  having  negative  values  of  the  displacements  shown,  but,  as  has  been 
discussed  earlier,  these  values  are  identical  to  the  extinction  efficiencies  of  systems  having 
positive  displacements. 

Figure  9  shows  the  extinction  efficiencies  as  a  function  of  the  inclusion  radius  when  the 
inclusion  is  at  the  edge  of  the  host  sphere  (aj  =  a2  d)  when  q  =  0®.  From  system  symmetiy" 
and  the  reciprocity  theorem,  Qextict)  =  Qext{—o:)  =  (5eit(180°  -  a).  A  comparison  of  Figures 
8  and  9  shows  that  there  is  a  greater  range  in  extinction  efficiencies,  and  therefore  a  greater 
sensitivity  of  the  efficiencies,  when  the  angle  of  incidence  is  varied  than  when  the  inclusion  is 
translated  in  the  direction  of  the  incident  beam. 

Examinations  of  the  above  figures  reveal  two  important  findings  for  the  nonconcentric  sphere 
system.  The  first  finding  is  that  the  extinction  efficiency  is  independent  of  the  sign  of  the 
displacement  of  the  inclusion  particle.  In  other  words,  the  system  gives  the  same  value  of  the 
extinction  efficiency  regardless  of  whether  the  inclusion  lies  in  the  forward  direction  (where 
the  fields  are  generally  stronger)  or  in  the  backscatter  direction  (where  the  fields  are  generally 
weaker).  We  also  see  that  the  closer  an  inclusion  is  to  the  surface  of  the  host  sphere,  the  greater 
the  effect  it  has  on  the  extinction  of  the  system.  This  last  point  is  veiy  important  in  terms  of 
resonance  suppression. 


5  CONCLUSION 


We  have  derived  solutions  for  the  total  field  when  an  incident  plane  wave  strikes  a  host 
sphere  containing  a  nonconcentric  spherical  inclusion.  Reduction  of  the  general  result  to  that 
of  the  concentric  sphere  and  the  Mie  sphere  verifies  the  validity  of  the  derivation.  We  have 
examined  the  scattering  as  a  function  of  a  number  of  system  parameters.  Of  particular  interest 
is  the  independence  of  the  extinction  efficiency  on  the  sign  of  the  displacement  of  the  inclusion 
particle.  We  also  see  that  the  closer  an  inclusion  is  to  the  surface  of  the  host  sphere,  the  greater 
the  effect  it  has  on  the  extinction  of  the  system. 
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Appendix:  NONCONCENTRIC  SPHERE  CODE 


Description  of  The  Program 


This  appendix  contains  the  FORTRAN  code  for  the  nonconcentric  spheres.  The  program 
consists  of  the  main  routine  and  twelve  subroutines.  Most  calls  to  the  subroutines  are  made 
from  the  main  program.  We  will  briefly  describe  each  subroutine  and  its  respective  variables. 
This  program  calculates  the  Mueller  scattering  matrix  elements  (5^),  and  the  extinction  (Qeif). 
scattering  {Qsca),  and  absorption  {Qahs)  efficiencies  for  a  set  of  initial  inputs.  The  input  data 
is  read  from  the  “input.in”  file.  All  the  variables  used  are  defined  in  a  separate  file  named 
“declare.def’.  The  reason  for  using  this  method  is  the  convenience  of  having  all  the  variables 
in  one  separate  file  for  comparison  and  debugging  purposes.  Explanation  of  some  important 
variables  are  given  in  the  subsections  below.  Also,  a  ^ort  descriptions  of  the  suteoutines  are 
presented. 

The  only  restriction  placed  on  this  program  is  that  we  cannot  allow  the  inclusion  radius 
to  eqilal  zero,  i.e.  02  7^  0.  The  reason  is  that  the  Neumann  functions  blow  up  if  the  argument 
is  zero.  If  the  user  wishes  to  have  no  inclusion,  the  user  can  let  the  index  of  refraction  of  the 
inclusion  equal  that  of  the  host  sphere,  i.e.  mj  =  m2. 


Subroutine  bessel  (n,np,bessj,x) 

This  subroutine  calculates  the  spherical  Bessel  function  of  the  first  kind  from  order  -1  to 
Older  n,  and  stores  the  results  in  the  array  bessj.  The  variable  n  is  related  to  the  size  parameter  x. 
If  X  is  large,  then  n  is  also  large.  The  variable  np  is  the  maximum  size  (the  physical  dimension) 
of  the  array.  Here  np  is  isize,  where  isize  is  defined  in  the  file  “deciare.def 

Input 

n:  The  largest  order  of  the  Bessel  function,  which  is  determined  in  the  main  routine  by 
the  definition  of  nbg. 

np:  The  array  size  defined  by  the  parameter  isize. 

x:  The  argument  of  the  Bessel  function. 

Output 

bessj:  This  double  precision  real  array  contains  the  values  of  the  spherical  Bessel  function 
of  the  first  kind  UniKai)] ,  from  order  -1  to  order  n. 
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Subroutine  C-bessel(n,np,bessj,x) 

This  subroutine  is  similar  to  the  Bessel  subroutine  except  that  it  takes  a  complex  argument 
X  instead  of  a  real  argument.  The  array  bessj  is  a  complex  double. 

Input 

n:  The  largest  order  of  the  Bessel  function,  which  is  determined  in  the  main  routine  by 
the  definition  of  nbg. 

np;  The  array  size  defined  by  the  parameter  isize. 
x:  The  complex  argument  of  the  Bessel  function. 

Output 

bessj:  This  double  complex  array  contains  the  complex  values  for  this  complex  spherical 
Bessel  function  of  the  first  kind  [jn(A:iai)j .  from  order  -1  to  order  n. 


Subroutine  newman(n^np,bessypc) 

This  subroutine  calculates  the  spherical  Bessel  function  of  the  second  kind  (or  the  Neumann 
function.  ni,{koai)),  from  order  -1  to  order  n  and  stores  the  results  in  the  double  precision  array 

bess}’. 

input 

n:  The  largest  order  of  the  Neumann  function,  which  is  determined  in  the  main  routine  by 
the  definition  of  nbg. 

np:  The  array  size  defined  by  the  parameter  isize. 
x:  The  argument  of  the  Neumann  function. 

Output 

bessy:  This  double  precision  real  array  contains  the  values  for  this  spherical  Neumann 
function  from  order  -1  to  order  n. 


Subroutine  C-newman(n,np,bessy,x) 

This  subroutine  is  similar  to  the  newman  subroutine,  except  that  it  takes  a  complex  argu¬ 
ment  X  instead  of  a  real  argument.  The  array  bessy  is  now  a  complex  double. 

Input 

n:  The  largest  order  of  the  Neuman  function,  which  is  determined  in  the  main  routine  by 
the  definition  of  nbg. 

np:  The  array  size  defined  by  the  parameter  isize. 
x:  The  complex  argument  of  the  Neumann  function. 

Output 

bessj:  This  double  complex  array  contains  the  complex  values  for  this  spherical  Neumann 
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function  i/7j,(A’]ai)j.  from  order  -1  to  order  n. 


Subroutine  riccati(nbg,hval,zeta,dzeta,x) 

This  subroutine  takes  as  input  the  Hankel  functions  of  the  first  or  second  kind  and  transform 
them  into  the  Riccati-Bessel  functions  and  their  derivatives.  The  transformation  follows  from 
the  definition  of  the  Riccati-Bessel  functions  shown  in  die  preceding  derivations. 

Input 

nbg:  The  largest  order  for  the  Hankel  and  Riccati-Bessel  functions, 
hval:  The  double  complex  array  which  contains  either  the  first  or  the  second 

[/i^(A:iai)]  Hankel  functions  as  defined  and  computed  in  the  main  routine.  The  values  stored  in 
this  array  range  from  order  -1  to  order  nbg. 

x;  The  double  complex  argument  of  the  spherical  Hankel  function. 

Output 

zeta:  This  is  the  Riccati-Bessel  function.  The  values  are  stored  from  order  zero  to  order 
nbg  into  this  double  complex  array. 

dzeta:  This  is  the  derivative  of  the  Riccati-Bessel  function,  the  results  are  stored  into  this 

auubic  i.uulpicx  diidv. 


Subroutine  nplgndr(m,nbr,iisize,legpol,x) 

Given  a  value  for  n%  where  m  <  nbr,  this  subroutine  calculates  the  normalized  associated 
Legendre  polynomials  [P”'(x)j  from  order  m  =  0  to  nbr  and  stores  the  values  in  the  array 
legpol. 

Input 

m:  The  integer  order  for  legpol.  which  ranges  from  0  to  nbg. 

nbr:  This  variable  has  a  maximum  value  of  nbg+1;  the  iteration  of  this  subroutine  goes 
from  m  to  nbr. 

iisize:  This  variable  is  defined  as  isize  +  1  (shown  in  the  “declare.def’  file). 

x:  This  argument  for  the  normalized  associated  Legendre  polynomials  is  the  cosine  of  the 
angle,  i.e.  x  =  cjos9. 

Output 

legpol:  This  is  the  one-dimensional  array  containing  the  calculated  values.  It  has  a  physical 
dimension  of  isize  +  1 . 


Subroutine  ludcmp(a,n,np,indx) 
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This  subroutine  decomposes  and  Tri-diagonalizes  a  matrix.  Given  a  matrix  a(l  :n,l  :n)  with 
physical  dimension  np,  where  np  is  defined  as  isiz2  in  the  “declare.def’  file,  this  subroutine 
replaces  it  by  the  LU  Decomposition  of  a  row-wise  permutation  of  itself.  The  procedure  is  an 
V?  process,  so  that  for  large  size  parameters,  the  computation  time  slows  down  dramatically.  If 
nbg  is  bigger  than  70  (which  corresponds  to  a  size  parameter  of  x  «  50),  one  should  change 
the  value  of  nmax  in  this  subroutine  so  that  one  always  has  nmax  >  2{nbg  +  1).  This  routine 
is  used  in  conjunction  with  the  subroutine  lubksb^  to  solve  linear  equations  or  invert  a  matrix. 
Note  that  we  have  modified  this  subroutine  to  accept  double  complex  arguments. 

Input 

a:  This  is  a  double  complex  n  x  n  matrix  given  by  equations  40  and  41  in  the  preceding 
derivation.  After  this  matrix  passes  through  this  subroutine,  it  will  be  destroyed  and  replaced 
with  its  tri-diagonal. 

n;  This  integer  value  is  equal  to  2(nbg+l).  The  parameter  nmax,  defined  within  this 
subroutine,  must  be  greater  than  or  equal  to  this  value. 

np:  This  is  the  physical  dimension  size  of  the  matrix  and  has  a  value  of  isiz2.  It  is  safe  to 
let  it  be  larger  than  n. 

Output 

a:  The  original  matrix  is  now  a  tri-diagonalized  matrix.  This  new  matrix  can  now  be  used 

in  ccnjunction  with  the  subroutine  lubksb  to  sob-e  the  linear  equation. 

indx:  This  is  an  array  which  records  the  row  permutation  effected  by  the  partial  pivoting. 


Subroutine  lubksb(aqi,np,indx,b) 

This  subroutine  solves  the  set  of  n  linear  equations  ax  =  b.  The  solution  vector  is  rep¬ 
resented  by  the  array  b.  and  the  solution  vector  x.  after  going  through  the  subroutine,  is  also 
stored  in  b.  The  variables  n,  np,  and  the  arrays  a  and  indx,  are  not  modified  by  this  routine  and 
can  be  left  in  place  for  successive  calls  with  different  right-hand  sides  b.  This  routine  takes  into 
account  the  possibility  that  b  will  begin  with  many  zero  elements,  so  it  is  efficient  for  use  in 
matrix  inversion.  We  have  also  modified  this  subroutine  to  accept  double  complex  arguments. 
Input 

a:  This  double  complex  matrix,  with  physical  dimension  np,  is  the  same  matrix  a  which 
has  undergone  the  tri-diagonalization  in  the  subroutine  ludcmp. 

n:  This  integer  value  is  equal  to  2(nbg+l).  The  variable  nmax.  defined  within  this  sub¬ 
routine,  must  be  greater  than  or  equal  to  this  value. 

np:  This  is  the  physical  dimoision  size  of  the  matrix  and  has  a  value  of  isiz2.  We  must 
keep  in  mind  that  np  must  be  greater  than  n. 

indx:  This  is  the  output  vector  which  records  the  row  permutation  effected  by  the  partial 
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pivoting  in  the  subroutine  ludcmp. 

b:  This  input  vector  is  the  right-hand  side  of  the  equation  ax  =  b,  i.e.  it  is  the  solution 
vector.  The  subroutine  will  solve  for  the  vector  x,  and  places  its  values  in  b. 

Output 

b:  The  vector  x  is  now  placed  in  this  array. 


Subroutine  mprove(a,alud,n,np,indx,b,x) 

This  subroutine,  which  has  been  modified  to  accept  double  complex  arguments,  improves 
a  solution  vector  x  of  the  linear  set  of  equations  ax  =  b.  It  assumes  that  the  vector  x  is  olf  from 
the  true  solution  by  bx,  and  it  recursively  solves  for  6x  and  improves  that  vector.  It  does  this 
by  comparing  the  two  matrices  a  and  alud,  with  the  solution  b  and  argument  x,  and  improves 
on  the  solution  of  x.  Note  that  inside  this  subroutine  is  the  parameter  nmax,  which  should  have 
the  same  value  as  that  given  in  the  subroutine  ludcmp. 

Input 

a;  This  nxn  matrix  is  the  original  matrix  before  it  has  undergone  any  changes  due  to  the 
result  of  the  subroutine  ludcmp.  All  the  values  are  in  their  original  state. 

aluJ.  Tois  n.  //.  matrix  is  the  one  that  underwent  the  tri-diagonalization  associated  v.ith 
the  ludcmp  subroutine.  It  is  the  output  matrix,  a,  given  in  the  ludcn^i  subroutine  above, 
n:  The  size  of  the  two  matrices,  and  the  size  of  the  vector  x,  b,  and  indx. 
np:  The  physical  dimension  of  the  two  matrices, 
indx;  The  ouqiut  vector  returned  by  tiie  subroutine  ludcmp. 
b:  The  solution  vector  of  the  equation  ax  =  b. 
x;  This  vector  is  the  output  from  the  subroutine  lubksb. 

Output 

x:  After  mprove,  x  has  been  modified  and  improved  to  a  new  set  of  values.  Note  that  this 
subroutine  can  be  called  as  often  as  one  wishes,  but  generally,  one  or  two  calls  are  enough. 


Function  find-g(nbg,nn,ang,cte,dte,ctm,dtm) 

This  function  calculates  the  asymmetry  parameter,  g.  and  returns  the  value  in  double 
precision  to  the  main  routine. 

Input 

nbg;  The  largest  order  of  the  Bessel  function  defined  in  the  main  routine, 
nn;  The  physical  size  of  the  arrays  used, 
ang:  The  incident  angle. 

cte:  This  is  the  array  which  contains  the  scattering  coefficients  for  the  TE  mode, 
dte:  This  too  contains  the  scattering  coefficients  for  the  TE  mode. 


A-5 


ctm:  This  array  contains  the  scattering  coefficients  for  the  TM  mode, 
dtm:  This  contains  the  scattering  coefficients  for  the  TM  mode. 
Output 

fmg-g;  Returns  as  a  double  precision  real  number. 


Subroutine  rotation! a,b,c,d,nbg,ang,cte,dte,ctm,dtm^,ij ) 

This  subroutine  is  used  in  conjunction  with  the  function  find-g  to  calculate  the  asymmetry 
parameter.  We  must  rotate  and  integrate  over  the  entire  sphere  in  order  to  find  the  asymmetry 
parameter. 

Input 

nbg:  The  largest  order  of  the  Bessel  function  used. 

ang:  The  incident  angle. 

nn:  The  physical  size  of  the  arrays  sent  from  the  main  routine. 

ij:  The  physical  size  of  the  arrays  sent  from  the  function  find-g. 

cte,  dte,  ctm,  dtm:  The  arrays  which  contain  the  scattering  coefficients  as  sent  in  from  the 
main  routine.  These  values  are  now  stored  into  the  arrays  a,  b,  c,  d  respectively. 

Output 

a,  b,  c,  d:  These  are  the  neu'  rotated  arrays  containing  the  scatttering  coefficients.  Now 
we  can  use  these  new  arrays  to  calculate  the  asymmetry  parameter. 


Function  factorial(n) 

This  function  calculates  the  factorial  of  a  given  integer.  Since  a  factorial  of  a  large  integer 
will  blow'  up.  we  define  this  fimction  to  take  the  sum  of  the  log  of  the  number.  This  will  ensure 
us  of  stability  in  our  calculations. 

Input 

n:  The  integer  value  sent  in  from  the  subroutine  rotation. 

Output 

factorial:  This  function  now  returns  a  double  precision  real  number. 


Declare,  def 

Within  the  main  routine  is  an  include  statement  which  refers  to  this  file.  The  parameter 
isize  is  presently  set  to  handle  a  size  parameter  of  up  to  50.  The  program  can  be  made  to  handle 
an  even  larger  size  parameter  by  increasing  isize  appropriately.  Note  that  if  the  size  parameter 
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of  the  run  is  greater  than  50,  one  must  change  isize  to  a  larger  number  to  include  all  memory 
allocations;  also,  the  parameter  nmax  in  the  subroutines  ludcmp  and  mprove  must  be  increased 
in  value.  The  other  parameter,  isiza,  is  presently  set  to  handle  up  to  360  degrees  by  increments 
of  9  degrees.  If  the  run  is  to  have  a  resolution  greater  than  1  degree  (which  can  be  changed  in 
the  input  file  with  the  variable  num),  isiza  must  be  changed  to  meet  that  specific  need. 
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c  . 

c  .  Scattering  program  for  an  inclusion  imbedded  inside  the 
c  .  host  droplet.  The  input  data  is  taken  from  the  file 

c  .  "input,  in"  and  the  variables  and  arrays  are  defined  in 

c  .  the  file  "declare.def. 

c  . 

implicit  double  precision  (a-h,o-z) 

INCLUDE  'declare.def 

open(unit^,file-input.in',status='old') 

rewind(4) 

open(unit=8, file-output. out') 

c  . 

c  .  read  in  the  input  values  from  "input. in"  . 

c  . 

read(4,*)  wavel 
read(4,*)  radl 
read(4,*)  rad2 
read(4,*)  distance 
readf4.*)  m_r1.m_i] 
read(4,*)  m_r2,m_i2 
read(4,*)  inc_ang 
read(4,*)  num 
read(4,*)  phi 

pi=1.0 

ccc=dcmplx(0.0,pi) 

pi=4.0*datan(pi) 

m_  1  =dcmpk(m_r  1  ,m_i  1 ) 

m_2=dcmpk(m_r2,m_i2) 

kO=2.*pLwavel 

kl=kO*m_l 

k2=k0*m_2 

xO=kO*radl 

xl=kl*radl 

x2=kl*rad2 

x3=k2*rad2 

c  nbg=  #  of  iterations  for  the  big  sphere  and  loops 

nbg=cdabs(xO)  +  4.*cdabs(x0)**.333  +  2.0 
print*,'the  value  of  nbg=',nbg 
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c  . 
c  . 


Calculate  Bessel  Functions  for  the  first  boundary  at  radl 
c  .  We  have  bessel  as  an  incident  field  of  the  sphere,  and  a 

c  .  hankel  as  the  scattered  field.  Their  arguments  are  xO  and  cxO 

c  .  Then  we  have  hankel  1  and  hankel2  inside  the  sphere;  and 

c  .  their  argument  is  kl *radl=xl . 

. . 

cxO=dreal(xO) 

call  bessel(nbg,isize,besj_o,cxO) 
call  newman(nbg,isize,temp,cxO) 
call  c_bessel(nbg,isize,hankel_o,xO) 
do  1 1  i=- 1  ,nbg 

hankel_o(i)=hankel_o(i)+ccc*temp(i) 

1 1  continue 

can  c_bessel(nbg,  isize, hankel  l_l,xl) 
can  c_bessel(nbg,isize,hankel2_l,xl) 
can  cnewman(nbg,isize,tempc,xl) 
do  12  i— l,nbg 

hankell_l(i)=hankell_l(i)  +  ccc*tempc(i) 
hankel2_l(i)=hankel2_l(i)  -  ccc*tempc(i) 

1 2  continue 


c  .  Now  we  match  boundary  at  rad2.  Here  we  have  two  hankels  in 
c  .  the  host  sphere,  and  their  argument  is  kl*rad2=x2 

c  .  Then  we  have  a  bessel  inside  the  inclusion  with  its  argument 

c  .  being  x3.  Note:  since  this  is  for  the  smaU  inclusion,  we 

c  .  must  take  care  to  have  spherical  bessel  fimctions  which  are 

c  .  within  the  Umits  of  that  size;  the  order  should  never  exceed 

c  .  the  size  of  the  sphere  because  it'U  blow  up  in  your  face 

. . 

nsm=cdabs(x2)  +  4.*cdabs(x2)**.333  +  2.0 
if(m_r2.gt.0.0)  cah  c_bessel(nsm,isize,besj_2,x3) 
can  c_bessel(nsm,isize,hankell_2,x2) 
can  c_bessel(nsm,isize,hankel2_2,x2) 
can  c_newman(nsm,isize,tempc,x2) 

do  13  i=-l,nsm 

hankell_2(i)=hankell_2(i)  +  ccc*tempc(i) 
hankel2_2(i)=hankel2_2(i)  -  ccc*tempc(i) 

13  continue 


Calculate  Riccati-Bessel  Functions  and  Derivatives  needed 
Note  that  zetal_l  means  the  1st  zeta  at  klal 


c  . 
c . 
c  . 


c  .  and  zeta2_l  means  the  2nd  zeta  at  klal 

c.  and  zeta  12  means  the  1st  zeta  at  kla2 

. . 

psi_o(- 1  )=cxO*besj_o(- 1 ) 
do  21  n=0.nbg 

psi_o(n)=cxO*besj_o(n) 
dpsi_o(n)=cxO*besj_o(n- 1 )  -n*besj_o(n) 

2 1  continue 

call  riccati(nbg,hankel_o,zeta_o,dzeta_o,xO) 
caD  riccati(nbg,hankell_l  ,zetal_l  ,dzetal_l  ,xl ) 
call  riccati(nbg,hankel2_  1  ,zeta2_l  ,dzeta2_  1  ,x  1 ) 
call  riccati(nsni,hankell_2,zetal_2,dzetal_2,x2) 
call  riccati(nsm,hankel2_2,zeta2_2,dzeta2_2,x2) 
call  riccati(nsm,besj_2,psi_2,dpsi_2,x3 ) 

. . 

c  .  Now  we  must  solve  for  the  Quality  factors  by  using  equations 
c  .  17  and  18;  These  Q's  relate  the  inclusion's  field  coeffs 

c  .  to  the  outter  sphere's  coefficients,  and  it  contains  infos  such 

c  .  as  the  inclusions's  radius  and  refractive  index 

c  .  Note  that  if  the  index  of  refraction  of  the  inclusion  is  nega- 
c  .  five,  then  we'll  take  that  to  mean  it  is  a  perfect  conductor, 

c  .  and  the  simplification  is  in  the  "else"  part  of  the  if-else 

c  .  statement.  We  cannot  use  index  >  100  because  it'll  blow  up 

. . 

if(m_r2.gt.0.0)  then 
do  25  n=0,nsm 

fact  1  =k  1  *dzeta2_2(n)*psi_2(n)-k2*zeta2_2(n)*dpsi_2(n) 
fact2=k2*zetal__2(n)*dpsi_2(n)-kl*dzetal_2(n)*psi_2(n) 
Q_r(n)=fact  1  /fact2 

fact3=k2*dzeta2_2(n)*psi_2(n)-kl*zeta2_2(n)*dpsi_2(n) 

fact4=kl*zetal_2(n)*dpsi_2(n)-k2*dzetal_2(n)*psi_2(n) 

Q_s(n)=fact3/fact4 

25  continue 

c  Now  to  fill  in  the  rest  of  the  array  with  a  constant  value.  This  is  a  must, 
do  26  n=nsm+l,nbg 
Q_r(n)=dcmplx(  1 .0,0.0) 

Q_s(n)=Q_r(n) 

26  continue 
else 

. . 

c  .  The  inclusion  is  a  perfect  conductor  . 

. . 

print*,The  inclusion  is  a  perfect  conductor' 
do  35  n=0,nsm 
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Q_r(n)=-zeta2_2(n)/zeta  1  _2(n) 
Q_s(n)=-dzeta2_2(n)/dzeta  1  _2(n) 

35  continue 

do  36  n=nsm+l,nbg 
Q_r(n)=dcmplx(  1 .0,0.0) 

Q_s(n)=Q_r(n) 

36  continue 
endif 

. . 

c  .  Calculate  translation  coefficients  when  the  separation  of  the 
c  .  center  to  center  radii  is  d.  We  will  first  find  the  bessel 
c.  function  which  describes  the  separation,  kid 

. . 

xd=(k0*m_l  )*distance 
if  (xd.eq.0.0)  then 
nsm=2 

else 

nsni=cdabs(xd)  +  4.*cdabs(xd)**.333  +2. 

endif 

call  c_bessel(3*nsm,isiz4,besj_d,xd) 


. . 

c  .  Begin  the  scalar  translation  coefficients  . 
c  .  keeping  in  mind  that  cuplo  =  c0,0(np) 
c  .  and  cuplo  1  =  c- 1 ,0(np)  for  starters 
. . 


do  110np=0,nbgM 
fiip=dfloat(np) 

cuplo(np)=besj_d(np)*dsqrt(2.  *fhp+ 1 .) 
cuplo  1  (np)=-cuplo(np) 

110  continue 

do  1 1 5  np=0,nbg*2 
cuplom(0,np)=cuplo(np) 

115  continue 

. . 

c  .  Now  for  some  acrobatics;  this  is  where  we  . 
c  .  perform  the  switcharoos  and  recursively  . 
c  .  solve  for  the  translation  coeff.  The 
c  .  result  is  stored  into  the  matrix  cuplom 
. . 


do  140n=l,nbg*2 
fh=dfloat(n) 
do  120  np=0,nbg*4-n 
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fiip=dfloat(np) 

c=cdsqrt(dcmplx((2.*fhp+l.)/(2.*fh-3.)))*(fo-l-) 

c=cuplol(np)*c 
cupio  1  (np)=cuplo(np) 
cuplo(np)=c 
120  continue 

do  1 30  np=0.nbg*4-n 
fiip=dfloat(np) 

c=-(lTip+ 1 .  )*dsqrt((2.  *fh- 1  .)/(2.  *fhp+3  .))*cuplo  1  (np+ 1 ) 
c=c+fiip*cdsqrt(dcmplx((2.*fh- 1 .  V(2.*fhp- 1  .)))*cuplo  1  (np- 1 ) 
c=c+cuplo(np) 

cuplo(np)=c*dsqrt((2.*fh+l.)/(2.*fhp+l.))/fri 
if  (np.le.nbg*2)  then 
cuplom(n,np)=cuplo(np) 
c  write(8,*)n,np,cuplom(n,np) 
endif 

1 30  continue 
140  continue 

c . 

c  .  Now  that  we  have  the  elements  in  the  matrix 
c  .  cuplom(n,np),  we  can  start  the  iterations  in 
c  .  m  to  get  the  rest  of  the  translation  coeff 
c  .  Then  we'll  have  the  complete  set  of  C(n,m)np 

c . 

do  299  m=0,nbg 

c  . 

c  .  But  first  we  need  to  calculate,  for  each  m,  the  . 

c  .  Legendre  Polynomials  for  incident  angle  theta  . 

c  . 

u=nbg+ 1 

x=dcos(pi*inc_ang/'  180.) 
mmm=m+l 

call  nplgndr(mmm,ii,isizel,legpoll,x) 
mmm=m 

call  nplgndr(mmm,u,isizel,legpol2,x) 
if  (m.gt.0)  then 
mmm=m-l 

call  nplgndr(mmni,ii,isizel,legpol,x) 
else 

do  165,iij— 2,isizel 
legpol(iii)=-legpoll  (iij) 

165  continue 
endif 
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do  1 80  i=m,nbg 
fi=dfloat(i) 

tau=dsqrt((fi+m+ 1  .)*(fi-ni))*legpoll  (i) 
pie=tau 

tau=(tau-dsqrt((fi+m)*(fi-m+ 1  .))*legpol(i))/2.0 
if  (dabs(x).gt.0.5)  then 

c  . . . 

c  .  The  angle  theta  is  between  0  and  60  degrees 

c  . 

pie=pie+dsqrt((fi+m)*(fi-m+l.))*legpol(i) 

pie=pie*.5/x 

else 

c  . 

c  .  theta  is  between  60  and  90;  use  this  so  the 

c  .  Polynomial  doesn't  blow  up  at  90  or  0  degrees 


pie=m*legpol2(i)/dsin(inc_ang*pi/ 1 80.) 
endif 

if  (i.gt.O)  then 

c  =  (ccc**i)/(fi*(fi+l.)) 

ba(i)=-c*pie 

aan)=c*tau 

endif 

1 80  continue 

. . 

c  .  That  ends  the  Calculation  for  the  Legendre  Polynomials 
c  .  for  each  value  of  m.  The  incident  field  coefficients 

c  .  for  the  TE  case  are  stored  in  aa(i)  and  ba(i);  the  TM 

c  .  case  will  be  done  later  by  flipping  the  TE  coefficients 

. . 

. . 

c  .  Now  to  find  the  rest  of  the  matrix  in  C(n,m)np.  For  each 
c  .  value  of  m  that  is  greater  than  0,  say  1,  we  can  find  the 
c  .  matrix  C(n,l)np  and  store  that  into  cuplom(n,np);  then 
c  .  the  next  iteration  of  m,  i.e.  m=2,  we  can  find  C(n,2)np 
c  .  by  already  having  C(n,l)np  stored  in  the  matrix  cuplom(n,np) 

. . 

if  (m.gt.0)  then 

do  155,  n=m-l,2*nbg-m+l 
do  154,  np=m-l,2*nbg-m+l 

cuplom  1  (n,np)=cuplom(n,np) 

1 54  continue 

155  continue 


do  160  n=m,2*nbg-m 
fh=dfloat(n) 
do  159  np=m.2*nbg-m 
fhp=dfIoat(np) 

c=dsqrt((fhp-m+l  .)*(ftip+m)*(2.*ftip+l  .))*cuplonil(n,np) 
cuplom(n,np)=c 

c=xd*dsqrt((fiip-m+2.)*(fiip-in+l.y(2.*fTip+3.)) 
cuplom(n,np)=cuplom(n,np)-c*cuplom  1  (n,np+ 1 ) 
c=xd*dsqrt((fhp+m)*(fiip+m- 1  .)/(2.*lhp-l .)) 
cuplom(n,np)=cuplom(n,np)-c*cuploml  (n,np- 1 ) 
c=dsqrt((fh-in+ 1  .)*(fo+ni)*(2.*fhp+l .)) 
cuplom(n,np)=cuplom(n,np)/ c 

159  continue 

1 60  continue 
endif 

c . 

c  .  Cuplom(n,np)  now  has  the  most  recent  values  . 
c  .  for  the  most  recent  value  of  m.  Note  that 
c  .  we  didn't  find  C(n,0)np  because  we  already 
c  .  have  those  values  stored  in  the  cuplom  matrix  . 
c . 

c . 

c  .  Now  that  we  have  the  matrix  which  contains  the  C(n,m)np  . 
c  .  we  can  go  ahead  and  find  the  vector  translation  coefficients  . 
c  .  A(n,m)np  and  B(n,m)np 

c . 

do  190  n=m,nbg 
pie=0.0 
tau=0.0 

do  195  np=m,nbg 
ftip=dfloat(np) 

tau=(fiip-ni+l  .0)*(ftip+m+l  .0) 
tau=tau/(  ( 2 .  *  fnp+ 1 . )  *  (2 .  *  fnp+3 . ) ) 
tau=dsqrt(tau) 

T  ransA=-xd/(fhp+ 1 .0) 

TransA=T  ransA*tau*cuplom(n,np+ 1 ) 

TransB=0.0 
if  (np.gt,0)  then 

pie=(fiip-ni)*(ftq>+ni) 
pie=pie/((2.*fiip- 1  .)*(2.*fhp+l .)) 
pie=dsqrt(pie) 

c=-xd*pie*cuplom(n,np- 1  )/fiip 
T  rans  A=T  ransA+c 

TransB=-ccc*xd*m*cuplom(n,np)/(fhp*(fhp+ 1 .)) 
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endif 

T  rans  A=cuplom(n,np)+T  ransA 
T  ransAm(  n,np)=T  ransA 
TransBm(n,np)=TransB 
195  continue 
190  continue 

numel=nbg-m+l 
numtot=numel+numel 
do  220  n=m,nbg 
i=n-m+l 

do  215  np=m,nbg 

. . 

c  .  The  Q's  have  the  index  of  the  inclusion  k2  inside  them 
c  .  and  the  TransAm  and  TransBm  have  the  displacement  d  inside 
c  .  them;  matrixlu  therefore  is  the  1  st  to  contain  both  pieces 
c  .  of  information. 

. . 

fact  1  =dzeta_o(n)*(zeta2_l  (n)+Q_r(np)*zetal_l  (n)) 
fact2=zeta_o(n)*(dzeta2_l  (n)+Q_r(np)*dzetal_l  (n)) 
fact3=dzeta_o(n)*(zeta2_l  (n)+Q_s(np)*zetal_l  (n)) 
fact4=zeta_o(  ni^f  dzeta2_U  n')+Q_s(np')*dzetal  _  1  (n)) 
ii=np-m+l 


c=TransAm(np,n)*(kO*factl  -kl*fact2) 

matrixlu(i,ii)=c 

matrix(i,ii)=c 


c=TransBm(np,n)*(kO*fact3  -  kl*fact4) 

matrixlu(i,ii+numel)=c 

matrix(i,ii+numel)=c 

c=TransBm(np,n)*(kl*factl  -  k0*fact2) 

matrixlu(i+numel,ii)=c 

matrix(i+numel,ii)=c 

c=TransAm(np,n)*(kl*fact3  -  k0*fact4) 
matrixlu(i+numel,ii+numel)=c 
matrix(i+numel,ii+numel)=c 
215  continue 
220  continue 


c.  THE  MATRIX  IS  FILLED!!  That  >vas  the  heart  of  the  problem  . 
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c . 

c  .  Now  we  can  load  up  the  solution  vectors  aa(n)  and  ba(n) 
c  .  and  start  the  LU  Decomposition  to  get  t(n)  and  u(n) 
c  .  Note  that  aa(n)  and  ba(n)  contain  the  incident  angle 

c . 

c . 

c  For  TE  case  . 

c . 

do  226  n=m,nbg 
i=n-m+l 

c=dzeta_o(n)*psi_o(n) 
c  =  kl *(c  -  dpsi_o(n)*zeta_o(n)) 
b(i)=aa(n)*c 
b(  i+numel)=ba(  n)  *  c 
bd(i)=b(i) 

bd(  i+numel)=b(  i+numel ) 

226  continue 
ii=isiz2 

call  ludcnip(matrixlu,numtot,ii,indx.dd) 
call  lubksb(matrixlu,numtot,ii,indx,b) 
call  mprove(niatrix,matrixlu,numtot,ii,indx,bd,b) 
do  227  n=m.nbg 
i=n-m+l 

TintTE(m,n)=b(i) 

UintTE(m,n)=b(i+numel) 

227  continue 

C  I  I  I  I  H-+  I  M  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  M  I  I  I  I  I  I  I  I  I  I 

c  .  Now  we  must  find  the  Scattering  coefficients  CscaTE  and  DscaTE, 
c  .  and  from  these  coeff  we  can  determine  the  scattering  matrices. 

. . 

do  250  n=m,nbg 
suml=0.0 
sum2=0.0 
do  260  np=m,nbg 

fact  1  ==dzeta2_  1  (n)+Q_r(np)*dzeta  1  _1  (n) 

fact  1  =T  ransAm(np,n)*T  intTE(m,np)*fact  1 

fact2=dzeta2_  1  (n)+Q_s(np)*  dzeta  1  _1  (n) 

fact2=TransBm(np,n)*UintTE(m,np)*fact2 

sum  1  =fact  1  +fact2+suml 

fact3=zeta2_  1  (n)+Q_r(np)*zeta  1  _1  (n) 

fact3=TransBm(np,n)*TintTE(m,np)*fact3 

fact4=zeta2_l(n)+Q_s(np)*zetal_l(n) 

fact4=TransAm(np,n)*UintTE(m,np)*fact4 

sum2=sum2+fact3+fact4 
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260  continue 

CscaTE(m,n)=(suml-aa(n)*dpsi_o(n))/dzeta_o(n) 

DscaTE(m,n)=(sum2-ba(n)*psi_o(n))/zeta_o(n) 

c  print*,n,CscaTE(ni,n),DscaTE(ni,n) 

250  continue 

c . 

cFOR  TM  CASE!!!  . 

c . 

do  228  n=m,nbg 
i=n-m+l 

c=dzeta_o(n)*psi_o(n) 

c  =kl*(c  -  dpsi_o(n)*zeta_o(n)) 
b(i)=ccc*ba(n)*c 
b(i+numel)=ccc*aa(n)*c 
bd(i)=b(i) 

bd(i+numel)=b(i+nuniel) 

228  continue 
ii=isi22 

call  lubksb(matrixlu,numtot,ii,indx,b) 

call  mprove(matrix,niatrixlu,niamtot,ii,indx,bd,b) 

do  229  n=m.nbg 
i=n-ni+l 

TintTM(m,n)=b(i) 

U  intTM(m,n)=b(i+numel) 

229  continue 

C  I  I  I  I  M  M  I  I  I  I  I  I  I  I  I  I  I  I  I  I  M  I  I  I  I  I  I  I  I  I  I  I  I  N  I  I  I  I  I  I  I  I  I  M  I  I  M' 

c  .  Now  we  must  find  the  Scattering  coelficients  CscaTM  and  DscaTM,  . 

c  .  and  Irom  these  coeflf  we  can  determine  the  scattering  matrices. 

. . 

do  25 1  n=m,nbg 
suml=0.0 
sum2=0.0 
do  261  np=m,nbg 

fact  l=dzeta2_l  (n)+Q_r(ip)*dzetal_l  (n) 

fact  1  =TransAm(np,n)*TmtTM(m,np)*fact  1 

fact2=dzeta2_l(n)+Q_s(np)*dzetal_l(n) 

fact2=TransBm(np,n)*UintTM(m,np)*fact2 

sum  1  =fact  1 +fact2+sum  1 

fact3=zeta2_l(n)+Q_r(np)*zetal_l(n) 

fact3=TransBm(np,n)*TintTM(m,np)*fact3 

fact4=zeta2_l(n)+Q_s(np)*zetal_l(n) 

fact4=TransAm(np,n)*UintTM(m,np)*fact4 

sum2=sum2+fact3+fact4 
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261  continue 

CscaTM(ni,n)=dpsi_o(n)*ccc*ba(n) 
CscaTM(m,n)={suml-CscaTM(ni,n))/dzeta_o(n) 
DscaTM(m,n)=ccc*aa(n)*psi_o(n) 
DscaTM(m,n)=(sum2-DscaTM(m,n))/zeta_o(n) 
c  print*, n,CscaTM(m,n),DscaTM(m,n) 

25 1  continue 


c  .  We  can  now  calculate  the  Efficiencies  (Qext,  Qsca,  Qabs) . 

c . 

do  510  n=l,nbg 

tau=CscaTE(m,n)*dconjg(CscaTE(m,n)) 
tau=tau+DscaTE(m,n)*dconjg(DscaTE(m,n)) 
pie=CscaTM(ni,n)*dconjg(CscaTM(m,n)) 
pie=pie+DscaTM(m,n)  *  dconjg(DscaTM(m,n) ) 
x=n*(n+ 1  .)*(tau+pie) 

s  1  =CscaTE(m,n)*dconjg(aa(n)) 
sl=sl+DscaTE(m,n)*dconjg(ba(n)) 
s2=CscaTM(ni,n)*dconjg(ccc*ba(n)) 
s2=s2+DscaTMfm.n')*dconjg('ccc*aa('n'0 
s3=n*(n+l.)*(s2+sl) 
if(m.eq.O)  then 
x=2.0*x 
s3=-2.0*s3 
else 

x=4.0*x 

s3=-4.0*s3 

endif 

qsca=qsca  +x 
qext=qext  +dreal(s3) 

5 1 0  continue 

c . . . 

c  .  That  finishes  one  loop  in  m,  now  we  move  to  the  next  m  and  repeat 
c  .  the  same  procedure  aU  over  again,  storing  our  internal  results 
c  .  into  T_te,  U_te,  T  tm,  U  tm,  and  storing  the  scattering  results 
c  .  into  C_te,  D_te,  C  tm,  D_tm. 

. . 

299  continue 

c . 

c  DONE ! !  The  sum  over  m  is  now  con^>lete! ! !  . 
c . 
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qext=qext/cdabs(xO)  *  *2 
qsca=qscaycdabs(xO)*  *2 
qabs=qext-qsca 

c  write(8,*)  real(distance),real(qext),real(qsca),real(qabs) 

print*,'Qext,  Qsca,  Qabs=',qext,qsca,qabs 

^  t*********************************************************** 

c  COMMENT  THIS  BLOCK  OUT  IF  YOU  DON'T  WANT  TO  FIND 

c  THE  ASYMMETRY  PAF;.AMETER.  THIS  FIND  G  SUBROUTINE 

c  IS  VERY  TIME  CONSUMING 
ii=isizel 
alpha=inc_ang 

c  gg=fmd_g(nbg,ii, alpha, CscaTE,DscaTE,CscaTM,DscaTM) 

gg=gg/xO*  *2/qsca 
print*,  'asym  parameter  =’,gg 

^  ilf;i:**t******************************************************* 


print*,  'Done  with  the  rough  stuff,  now  to  find  the  Mueller  elements' 
print*  '************************************************' 


. . 

c  .  Here's  the  meat  of  the  problem.  We  will  calculate  the  an^litude 
c  .  scattering  matrices  for  the  fields  that  we've  found.  The  angle 
c  .  theta  is  our  observation  angle.  Num  is  the  number  of  angles  we 
c  .  want  in  determining  the  scattering  intensity.  The  rest  of  this 
c  .  should  look  almost  like  the  Mie  Code  found  in  Bohren  and  Hufiman 
c  .  We  have  here  3  nested  loops.  The  first  one  is  to  determine  the 
c  .  angles  we  want  to  look  at,  the  second  one  is  the  sum  over  m,  and 
c  .  the  third  one  sums  over  n. 

c  .  If  we  want  to  have  the  observation  angle's  resolution  of  1  degree 
c  .  then  set  num=360,  if  we  want  finer  resolution,  let  num  be  a  bigger 
c  .  number  such  as  720  (to  have  resolution  of  0.5degrees) 

. . 


phi_o=df]oat(phi*pi/ 1 80.0) 
i=0 

dang=360.0/num 

b_angle=inc_ang+ 180. 
f_angle=inc_ang 
do  399  j=l,num+l 
angle=d_ang*(j- 1 ) 
angledd=angle 

if((angle.le.  1 80.0).and.(angle.ge.0.0))  then 
print*,'Obs  ANGLE  IS  ',angle 
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if^angle.gt.  180.0)  then 
angle=3 60 . 0-angle 
phi=pi  +  phi_o 
else 

phi=phi_o 

endif 

i=i+l 

slTE(i)=0. 

s3TE(i)=0. 

s2TM(i)=0. 

s4TM(i)=0. 

theta=angle*pi/ 1 80. 
x=dcos(  theta) 

do  330  m=0,nbg 

.  The  next  few  lines  will  calc  the  Legendre  polynomial. 
.  at  the  angle  theta 


ii=nbg+l 

mmm=m 

call  nplgndr(mmni,ii,isizel,legpol2,x) 
nimm=m+l 

caU  nplgndr(mmni,ii,isize  1  ,legpol  1  ^x) 
if  (m.eq.O)  then 

do  305ij=-2,isizel 

legpol(ij)=-legpoll(ij) 

continue 

else 

mmm=m-l 

caU  nplgndr(mmni,ii,isizel,legpol,x) 
endif 

do  310  n=m,nbg 


.  Calculate  the  corresponding  pie  and  tau  . 


tau=sqrt((n+m+ 1  .)*(n-m))*legpoll  (n) 
pie=tau 

tau— (tau-sqrt((n+m)*(n-m+l.))*legpol(n))/2. 
if  (dabs(x).gt.0.5)  then 
pie=pie+sqrt((n+m)*(n-m+ 1  .))*legpol(n) 
pie=pie*.5/x 
else 


pie=m*legpol2(n)/dsin(  theta) 
endif 

. . 

c  .Done  with  finding  Tau  and  Pie! 
c  .Now  to  find  the  scattering  amplitude  Sl,S2,S3,and  S4  . 
. . 


ifact==(-ccc)**n 
if(m.eq.O)  then 

fact  1  =DscaTE(0,n)*pie+CscaTE(0,n)*tau 
fact2=CscaTM(0,n)*pie+DscaTM(0,n)*tau 
fact3=CscaTE(0,n)*pie+DscaTE(0,n)*tau 
fact4=DscaTM(0,n)*pie+CscaTM(0,n)*tau 
sl=ifact*factl 
s2=ifact*fact2 
s3=ifact*fact3 
s4=ifact*fact4 
else 

fact  1  =DscaTE(in,n)*pie+CscaTE(ni,n)*tau 
fact  1  =fact  1  *2.0*dcos(m*phi) 
fact2=CscaTM(m,n)  *pie+DscaTM(in,n)  *tau 
fact2=fact2  *2.0*dcos(m*phi) 
fact3=CscaTE(m,n)*pie+DscaTE(m,n)*iau 
fact3=fact3  *2.0*ccc*dsin(m*phi) 
fact4=DscaTM(m,n)*pie+CscaTM(m,n)*tau 
fact4=fact4*2.0*ccc*dsin(m*phi) 
sl=ifact*factl 
s2=ifact*fact2 
s3=ifact*fact3 
s4=ifact*fact4 
endif 

slTE(i)=slTE(i)+sl 
s2TM(i)=s2TM(i)+s2 
s3TE(i)=s3TE(i)+s3 
s4TM(i)=s4TM(i)+s4 
3 1 0  continue 

. . 

c  .  Finished  summing  over  n,  now  need  to 
c  .  sum  over  m  in  order  to  have  the  complete  . 
c  .  solution  for  the  Scat  Amp  Matrix 

. . 

330  continue 


slTE(i)=slTE(i) 

s2TM(i)=-ccc*s2TM(i) 
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s3TE(i)=-ccc*s3TE(i) 

s4TM(i)=s4TM(i) 

c . 

c  .  Done  summing  over  m.  Now  we  have  one 
c  .  complete  set  of  scat  amp  S(i) 
c . 

c  -I-++I  I  M  I  M  I  I  M  I  I  I  I  I  M  I  I  M  I  I  M  f-l  I  I  I  M  I  f  I  I  I  I  II  I  I  I  I  I  I  I  I  M  I  I  I  I  I  M  I  I 

c  .  Lastly,  we  can  now  use  the  scattering  amplitude  matrix 

c  .  and  solve  for  the  Mueller  Matrix  elements.  Once  we  have 
c  .  the  MM.  we  have  all  the  informations  we  need. 

C  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  4  ■f+ 1  I  I  I  I  I  I  +4-l-++-H-+-H-t-+ 1  I  )  M  +H  I  I  I  I  I  M  I  I  I  I  I  I  ■^^++4 

do  410  jj=l,4 
mm(l,jy)=0.0 
mm(2,]j)=0.0 
mm(3,jj)=0.0 
mm(4,jy)=0.0 
410  continue 

fact  l=s  1  TE(i)*dconjg(s  1  TE(i))/2. 
fact2=s2TM(i)*dconjg(s2TM(i))/2. 
fact3=s3TE(i)*dconjg(s3TE(i))/2. 
fact4=s4TM(i)*dconjg(s4TM(i))/2. 
c=s2TM(i)*dconjgfs3TE(i'» 
cc=slTE(i)*dconjg(s4TM(i)) 
mm(  1 , 1  )=dreal(fact  1  +fact2+fact3+fact4) 
mm(  1 ,2)=dreal(fact2-fact  1  -fact3+fact4)/mm(  1,1) 
mm(  1 ,3  )=dreal(c  +  cc)/mm(  1,1) 
mm(  1 ,4)=dimag(c  -cc)/mm(  1,1) 

mm(2, 1  )=dreal(fact2-fact  1  +fact3-fact4 ) 
mm(2,2)=dreal(fact2+fact  1  -fact3-fact4)/mm(  1,1) 
mm(2,3)=dreal(c  -cc)/mm(l,l) 
mm(2,4)=dimag(c+cc)/mm(  1,1) 

fact  1  =s2TM(i)*dconjg(s4TM(i)) 

fact2=slTE(i)*dconjg(s3TE(i)) 

fact3=sl  TE(i)*dconjg(s2TM(i)) 

fact4=s3TE(i)*dconjg(s4TM(i)) 

c=s2TM(i)*dconjg(s  1  TE(i)) 

cc=s4TM(i)*dconjg(s3TE(i)) 

mm(3, 1  )=dreal(fact  1  +fact2) 

mm(  3 ,2)=dreal(fact  1  -fact2)/mm(  1,1) 

mm(3,3)=dreal(fact3+fact4)/mm(  1,1) 

mm(3,4)=dimag(c  +  cc)/mm(  1,1) 

fact  1  =s4TM(i)*dconjg(s2TM(i)) 
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fact3=slTE(i)*dconjg(s2TM(i)) 
mm(4, 1  )=dmiag(fact  1  +fact2) 
nim(4,2)=dimag(fact  1  -fact2)/mm(  1,1) 
mm(4,3)=dimag(fact3-fact4)/imn(  1,1) 
nim(4,4)=dreal(fact3-fact4)/mm(  1,1) 
print*,real(angledd),real(mm(  1,1)) 
write(8,*)int(angle),real(mm(  1 , 1  )),real(nim(  1 ,2)), 

+  real(nini(3,3)),real(nim(3,4)) 
endif 

c - 

c  This  endif  statement  lets  us  select  the  angles  we  want  to  look  at  . 

c - 

399  continue 

c  . 

c  .  Done  with  calculating  the  amplitude  coefficients  . 

c  . 


print*.  That  is  all.  Folks! ' 

999  close(4) 
close(8) 
stop 
end 
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This  is  the  Input  File 


0.6328 

The  free  space  wavelength  (wavel) 

0.525000 

The  host  radius  in  microns  (radl) 

0.125 

radius  of  the  inclusion  (rad2  ) 

0.1000 

center-center  separation  distance  (note:  d+rad2  <  radl) 

1.550  O.Od-1 

refractive  index  of  host  (ml ) 

1.550  O.Od-1 

refractive  index  of  inclusion  (m2) 

0.0 

incident  angle  of  beam  (inc  ang) 

20 

number  of  obs.  angles  between  0  and  360  (num) 

0.0 

angle  phi  (relative  to  plane) 
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c  . 

c  DECLARE.DEF 

c 

c  .  This  is  the  declaration  for  the  variables  which  will  be 

c  .  used  in  the  scattering  program.  All  the  arrays  and  other 

c  .  variables  can  be  declared  and  changed  in  here, 
c  . 

c  The  parameter  ♦isize*  is  the  physical  dimension  of  the  arrays.  We 
c  should  make  it  at  least  as  big  as  nbg  (suggest  nbg+2).  Also, 
c  you  should  make  sure  that  the  parameter  nmax  in  the  subroutines 
c  mprove  and  ludcmp  should  be  at  least  2*isize.  If  you  don't  care 
c  about  the  scattering  matrix  elements,  set  isiza=l ;  otherwise,  set 
c  isiza=361  (one  for  each  degree  from  0  to  360  degrees). 

parameter  (isize=200,  isiza=361) 
parameter  (isize  1  =isize+ 1 ) 
parameter  (isiz2=2 ♦isize,  isiz4=4*isize) 


c  For  the  MEDIUM,  which  in  our  case  is  just  plain  old  air. 
double  precision  besj_o(-l  :isize),temp(-l  :isize) 
double  precision  psi_o(-l  :isize).dpsi_o(-l  risizei 
complex  ♦  1 6  hankel_o(-l  :isize),tempc(-l  risize) 
complex  *  1 6  zeta_o(-l  :isize),dzeta_o(-l  :isize) 

c  For  the  HOST  SPHERE,  all  the  arrays  need  be  at  least  nbg 
complex  *16  hankell_l(-l:isize) 
complex  *16  zetal_l(-l:isize),dzetal_l(-l:isize) 
complex  *  1 6  hankel2_l  (-1  :isize) 
complex  ♦  1 6  zeta2_l  (-1  :isize),dzeta2_l  (-1  :isize) 

c  For  the  INCLUSION  SPHERE 

coit^jlex  *  1 6  besj_2(-l  ;isize) 

complex  *16  psi_2(-l  :isize),dpsi_2(-l  :isize) 

complex  *16  hankell_2(-l  :isize) 

complex  *16  zetal_2(-l:isize),(lzetal_2(-l  risize) 

complex  *  1 6  hankel2_2(-l  risize) 

complex  *16  zeta2_2(-lrisize),dzeta2_2(-l  risize) 

c  For  the  separation,  and  besj_d  need  be  at  least  4*nbg 
complex  *16  xd,  besj_d(-l  risiz4) 

c  For  the  refractive  indices  and  k-vectors 

double  precision  m_rl  ,m_r2,m_il,m_i2 
complex  *  1 6  k0,kl  ,k2,x0,xl  ,x2,x3,m_l  ,m_2 
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c  For  the  Internal  Coefficients 

complex  *16  TmtTE(0;isizel,0:isizel),  UintTE(0:isizel,0:isizel) 
con^lex  *16  TintTM(0:isizel,0:isizel),  UintTM(0:isizel,0:isizel) 

c  For  the  SCATTERING 

conqjlex  *16  CscaTE(0:isizel,0:isizel),DscaTE(0:isizel,0:isizel) 
complex  *16  CscaTM(0:isizel,0:isizel),DscaTM(0;isizel,0:isizel) 
conplex  *  1 6  Q_r(0;isize),Q_s(0:isize) 
complex  *16  slTE(isiza),s2TM(isiza),s3TE(isiza),s4TM(isiza) 
complex  *16  sl,s2,s3,s4,suml,sum2,factl,fact2,fact3,fact4 

c  For  the  NORMALIZED  ASSOCIATED  LEGENDRE  POLYNOMIALS 
double  precision  legpol(-2;isizel  ),legpoll  (-2:isizel ) 
double  precision  legpol2(-2;isizel),  tau,pie 
double  precision  inc_ang,x 

c  For  the  TRANSLATION  COEFFICIENTS 

complex  *  1 6  cuplo(-l  :isiz4),  cuplo  1  (-1  :isiz4) 

complex  *16  cuplom(0:isiz2,0:isiz2),  cuploml(0:isiz2,0:isiz2) 

complex  *16  TransA,TransB,TransAm(0:isize,0:isize) 

complex  *  1 6  TransBm(0:isize,0:isi2e) 

complex  *  1 6  c.cc.ccc.ifact 

c  For  the  MATRICES  IN  THE  LU_DECOMPOSTION 
dimension  indx(isiz2) 

complex  *  1 6  matrix(isiz2,isk2),  matrixlu(isiz2,isiz2) 
cort^lex  *16  b(isiz2),bd(isiz2),aa(0:isize),ba(0:isize) 

c  MUELLER  MATRIX 

double  precision  mm(4.4) 
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c 

c 


:Jc**3ic:t:*:*:*is**3*:**5|t**  *♦♦♦*****♦**♦****♦♦**  *******  *********** 

**********  these  are  the  subroutines  **************** 
******************************************************** 

subroutine  bessel(n,np,bessj,x) 


c  .subroutine  to  calculate  spherical  bessel  functions  of  the  first 
c  .kind  from  order  -1  to  order  n,  the  result  is  stored  in  the  array 
c  .bessj.  n  is  related  to  the  size  parameter  (nbg),  and  np  is  the 
c  .physical  dimension  of  the  bessj  array  (isize).  x  is  a  real  argument 


implicit  double  precision  (a-h,o-z) 
parameter(iacc=500,bigno= 1  d  1 00,bigni=  1  d- 1 00) 
integer  n,np 

double  precision  bessj(- 1  :np),x,xx,pi 
double  precision  bjp,bj,factor,ax 

pi=l  .0 

pi=4.*datan(pi) 
xx=>:-2*pi*intrx^27pil 
if  (n.lt.3)  then 
n=3 
endif 

. . . . 

c  .  If  X  is  less  than  ru  we  must  use  the  downward  recurrence  relation, 
c  .  otherwise  our  results  for  bessj  will  be  unstable.  In  fact,  we  will 
c  .  use  the  downward  alot  more  than  the  upward  because  our  x  will 
c  .  often  be  less  than  n,  i.e.  the  size  parameter  is  less  than  the  order 


if  (x.lt.n)  then 
ax=dabs(x) 

axx=x-2*pi*int(ax/2./pi) 
if  (ax.eq.O)  then 
bessj(0)=l. 
do  10  i=l,  n 
bessj(i)=0. 

10  continue 
else 

m=2*((n+int(sqrt(float(iacc*n))))/2) 


c  .  m  is  the  highest  order  to  calculate;  the  higher  m  is  . 
c  .  the  more  accuracy  we'll  get,  but  we'll  keep  n  terms  . 
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do  15j=0,n 
bessj(j)=0. 

1 5  continue 
bjp=0. 

bj=l. 

do  20 

bjm=(2.  *j+ 1  .)/ax*bj-bjp 

bjp=bj 

bj=bjm 

if(dabs(bj).gt.bigno)  then 

. . 

c  .  If  bj  is  huge,  we  must  renormalize  to  prevent  overflow  . 

. . 

bj=bj*bigni 
bjp=bjp*bigni 
do  17  i=0.n 

bessj(i)=bessj(i)*bigni 
1 7  continue 

endif 

if  (j.le.n)  bessj(j)=bjp 
20  continue 
bcssj(0)=bj 
sum=0. 
do  25  j=0,n 

suin=sum+(2.*j+ 1  .)*bessj(j)*bessj(j) 

25  continue 
do  30  j=0,n 

bessj(j)=bessj(j)/dsqrt(  sum) 
if  (x.lt.0.and.mod(n,2).eq.l)  bessj0')=-bessj0‘) 
30  continue 

bessj(- 1  )=dcos(xx)/x 
endif 
else 

. . 

c  .  Use  upward  recurrence  for  when  x  is  greater  than  n  . 


bessj(- 1  )=dcos(xx)/x 
bessj(0)=dsin(xx)/x 
do  40  i=0,n-l 

bessjXi+ 1  )=bessj(i)*(2.*i+l  .)/x-bessj(i- 1 ) 
40  continue 
endif 
return 
end 
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^  **************************************************** 
subroutine  c  bessel(n,np,bessj,x) 

. . 

c  .subroutine  to  calculate  conplex  spherical  bessel  functions  of  the  . 
c  .first  kind  from  order  -1  to  order  n.  The  result  is  stored  in  array  . 
c.bessj.  This  is  just  like  the  bessel  function,  except  we  now  have  . 
c  .complex  arguments,  i.e.  x  is  a  complex  variable 
. . 


parameter(iacc=500,bigno= 1  d  1 00,bigni- 1  d- 1 00) 

inplicit  double  precision  (a-h,o-z) 
double  precision  pi,xx 
integer  n,np 

complex  *16  x.xxxx,sum 
complex  *  16  bessj(-l  :np) 
complex  *16  bjp,bj,factor,bjm 

if  (n.lt.3)  then 
n=3 
endif 

xx=cdabs(x) 

pi=1.0 

pi=4.*datan(pi) 

xxxx=x 

if  (xx.gt. 32760)  then 

xxxx=x-2*pi*int(x/2./pi) 

endif 

if  (xx.lt.n)  then 

c . 

c  .  Use  downward  recurrence  . 

c . 

if  (xx.eq.O)  then 
bessj(0)=l. 
do  10  i=l,n 
bessj(i)=0. 

10  continue 

else 

m=2*((n+int(sqrt(float(iacc*n))))/2) 
do  1 5  j=0,n 
bessj(j)=0. 

1 5  continue 

bjp=dcmpk(0.d0,0.) 
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bj=dcmplx(l.dO,l.) 
do  20 

bjm={2.*j+l.)/x*bj-bjp 

bjp=bj 

bj=bjm 

xxx=cdabs(bj) 
if(xxx.gt.bigno)  then 

c . 

c  .  Renormalize  to  prevent  overflow  . 

c . 

bj=bj*bigni 
bjp=bjp*bigni 
do  1 7  i=0,n 

bessj  ( i)=bessj(i)  *bigni 
1 7  continue 

endif 

if  Ci-le.n)  bessj(j)=bjp 
20  continue 
bessj(0)=bj 
sum=0. 
do  25  j=0,n 

sum=sum+('  2 .  *j+ 1 .  i  *bessj(i')*bessj(i') 
25  continue 
do  30  j=0,n 

bessj(j)=bessj(j)/cdsqrt(sum) 

30  continue 

bessj(- 1  )=cdcos(xxxx)/x 
endif 

else 

c . 

c  .  Use  upward  recurrence  . 

c . 

bessj(- 1  )=cdcos(xxxx)/x 
bessj(0)=cdsin(xxxx)/x 
do  40  i=0.n- 1 

bessj(i+ 1  )=bessj(i)*(2.*i+l  )/x-bessj(i- 1 ) 
40  continue 

endif 
return 
end 


^  ***************************=*=************************ 
subroutine  newman(n,np,bessy,x) 

. . 

c  .  subroutine  to  calculate  the  spherical  bessel  function  of  the 
c  .  second  kind  (or  Newmann  function)  from  order  -1 
c  .  to  order  n;  result  stored  in  array  bessy.  bess_y  is  the  common  . 
c  .  notation  for  the  Newmann  function. 

. . 


implicit  double  precision  (a-h,o-z) 
integer  n,np 

double  precision  x,xx,pi 
double  precision  bessy(-l  :np) 

if  (n.lt.3)  then 
n=3 
endif 
pi=1.0 

pi=4.*datan(pi) 

xx=x-2.*pi*int(x/2./pi) 

. . . . •— 

c  .Must  use  this  xx  because  it's  a  value  between  0  and  2pi.  and  it . 

c  .works  as  the  argument  for  the  sin  and  cos;  if  we  just  stuck  in  . 

c  .X  instead,  it  won't  work. 

. . 

bessy(- 1  )=dsin(xx)/x 

bessy(0)=-dcos(xx)/x 

bessy(  1  )=-dcos(xx)/x/x-dsin(xx)/x 

do  10  i=2,n 

bessy(i)=(2.*(i-l .)+!  •yx*bessy(i-l)-bessy(i-2) 

1 0  continue 

return 

end 

^  ************************************************** 

subroutine  c_newman(n,np,bessy,x) 

. . . . 

c  .  subroutine  to  calculate  the  complex  spherical  bessel  function 

c  .  from  order  0  to  order  n,  the  result  is  stored  in  array  bessy. 

c  .  This  is  the  Newman  function,  aka,  bessel  function  of  order  2. 

. . 


implicit  double  pfecision(a-h,o-z) 
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integer  n,np 

complex  *16  x,xx,bessy(-l:np) 
double  precision  pi 

if  (n.lt.3)  then 
n=3 
endif 
pi=1.0 

pi=4.*datan(pi) 

xx=x-2.*pi*int(x/2./pi) 

bessy(- 1  )=cdsin(xx)/x 

bessy(0)=-cdcos(xx)/x 

bessy(  1  )=-cdcos(xx)/x/x-cdsin(xx)/x 

do  1 0  i=2,n 

bessy(i)=(2.*(i-l  .)+l  .yx*bessy(i-l)-bessy(i-2) 
1 0  continue 
return 
end 


subroutine  npigndrt m.nbr.iisize.legpol.x) 

. . 

c  .  Calculate  normalized  ass.  legendre  pol.  from  1=0  to  nbr  using  value  of 
c  .  -l<=x<=l  Note  that  everything  here  is  double  precision  because  other  . 
c  .  wise  weVe  got  a  runaway  solution. 

. . 


implicit  double  precision  (a-h,o-z) 
double  precision  x 
integer  m,nbr,iisize 
double  precision  legpol(-2:iisize) 
double  precision  pmm,somx2,fact,pll,pmmpl 
mm=mod(m,2)*(-2)+l 
m=abs(m) 
do  5  i=-2,m 
legpol(i)=0. 

5  continue 
pmm=  1. 
if  (m.gt.0)  then 

somx2=sqrt((  1  .-x)*(  1  .+x)) 

fact=l. 

do  1 1  i-l,m 

pmm— pmm*fact*somx2/sqrt(real(i)) 
fact=fact+2. 
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1 1  continue 

do  110  i=m+l,2*m 

pmm=pmm/sqrt(real(i)) 

1 1 0  continue 

pmm=pmm*  sqrt(2 .  *  m+ 1 . ) 
endif 

legpol(m)=pmm*nim 
if  (nbr.gt.m)  then 

pmmp  1  =x*sqrt(2.  *m+3  .)*pmm 
legpol(m+ 1  )=pmmp  1  *nim 
endif 

if  (nbr.gt.m+1)  then 
do  12  ll=m+2,nbr 

pll=x*sqrt(2.*ll- 1  .)*pmmp  1 
pll=pU-sqrt((U+m-l.)*(U-m-l.y(2.*U-3.))*pmm 

pl]=pll/sqrt((U-m)*(U+m)/(2.*Il+l.)) 

pnim=pnimpl 

pmmpl=pll 

legpol(U)=pll*mni 

12  continue 
endif 
m=m*mm 
return 

end 


subroutine  riccati(nbg,hval,zeta,dzeta,x) 
integer  nbg,i 

complex  *16  hval(-l:nbg),zeta(-l:nbg),dzeta(-l:nbg),x 

zeta(- 1  )=x*hval(- 1 ) 
do  10  i=0,nbg 

zeta(i)=x*hval(i) 
dzeta(i)=x*hval(i-l)  -i*hval(i) 

10  continue 
return 
end 
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subroutine  ludcmp(a,n,np,indx,d) 

c . 

c  .  given  nxn  matrix  a  with  physical  dim  NP,  this  routine  replaces  it 
c  .  by  the  LU  decomposition  of  a  rowwise  permutation  of  itself.  A  and 
c  .  n  are  input.  A  is  output.  INDX  is  an  output  vector  which  records  the 
c  .  row  permutation  effected  by  the  partial  pivoting;  D  is  output 
c  .  as  +-1  depending  on  whether  the  number  of  row  interchanges 
c  .  was  even  or  odd.  This  routine  is  used  in  combination  with  Lubksb 
c  .  to  solve  linear  equations  or  invert  a  matrix. 

. . 


complex  *  1 6  tiny 
parameter  (nmax=150) 

c . 

c  .  this  nmax  value  should  be  at  . 
c  .  least  twice  the  value  of  nbg 

c . 

complex  *16  a(np,np),suni,cdum,vv(nmax),aamax 
dimension  indx(n) 

tin^-dcmplxf  1  .Od-80.1  .Od-80’) 

d=l. 

do  12  i=l,n 

aamax=dcmplx(0.d0,0.) 
do  1 1  j=l,n 

if  (cdabs(a(ij)).gt.cdabs(aamax))  aamax=a(i,j) 

1 1  continue 

if  (cdabs(aamax).eq.O)  pause  'singular  matrix* 
vv(i)=l./aamax 

1 2  continue 
do  19j=l,n 

do  14i=l,j-l 
sum=a(i,j) 
do  13k=l,i-l 

sum=sum-a(  i,k)  *a(k,j) 

1 3  continue 
a(ij)=sum 

14  continue 
aamax=dcnpk(0.d0,0.) 
do  1 6  i=j,n 

sum=a(i,j) 
do  15  k=l,j-l 

sum=sum-a(i,k)*a(k,j) 

1 5  continue 
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a(i,j)=sum 

cdum=w(i)*sum 

if  (cdabs(c(ium).ge.cdabs(aamax))  then 
imax=i 
aamax=cdum 
endif 

16  continue 

if  (j.ne.imax)  then 
do  17k=l,n 

cdum=a(iniax,k) 

a(iniax,k)=a(_j,k) 

a(j,k)=cdum 

1 7  continue 
d=-d 

w(inriax)=w(j) 

endif 

indx(j)=imax 

if(cdabs(a(j,j))-eq.O)  a(j,j)=tiny 
if  (j.ne.n)  then 
cduni=l./a(j,j) 
do  18  i=j+l,n 

1 8  continue 
endif 

19  continue 
return 
end 

^  ^i^i^idt:***^****************************************** 

subroutine  lubksb(a,n,np,indx,b) 

. . . . 

c  .  Solves  the  set  of  n  linear  equations  A.X=B.  Here  A  is  input,  not 
c  .  as  the  matrix  A  but  rather  as  its  LU  decomp,  from  LUDCMP.  INDX 
c  .  is  input  as  the  permutation  vector  returned  by  LUDCMP.  B  is  input 
c  .  as  the  right-hand  side  vector  B,  and  returns  with  the  solution 
c  .  Vector  X.  A,N,NP  and  INDX  are  not  modified  by  this  routine  and 
c  .  can  be  left  in  place  for  successive  calls  with  different  right-hand 
c  .  sides  B.  This  routine  takes  into  account  the  possibility  that  B 
c  .  can  begin  with  many  zero  elements,  so  it  is  efficient  for  use  in 
c .  matrix  inversion. 

. . 

complex  *16  a(np,np),b(n),sum 
dimension  indx(n) 
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ii=0 

do  12  i=l,n 
ll=indx(i) 
siim=b(ll) 
b(U)=b(i) 
if  (ii.ne.O)  then 
do  1 1 

sum=sum-a(i,j)*b(j) 

1 1  continue 

else  if  (cdabs(suni).ne.O.)  then 
ii=i 
endif 
b(i)=sum 

12  continue 

do  14i-n.l,-l 
sum=b(i) 
if  (i.lt.n)  then 
do  13j=i+l,n 

sum=sum-a(i,j)*b(j) 

13  continue 
endif 

bfi)=sum/afi.i') 

14  continue 
return 
end 


************************************************** 


subroutine  mprove(a,alud,n,np,indx,b,x) 

. . 

c  .  improves  a  solution  vector  X  of  the  linear  set  of  equations  A*X=B  . 
c  .  The  matrix  A,  and  the  vectors  B  and  X  are  input,  as  is  the  dim 
c  .  n.  Also  input  is  alud  the  ludcomp  of  A  as  returned  by  ludcmp,  and  . 
c  .  the  vector  indx  also  returned  by  that  routine.  On  output,  only 
c  .  X  is  modified,  to  an  improved  set  of  values. 

. . 


parameter  (nmax=150) 

. . 

c  .  all  the  values  here  should  be  the  same  as  those 
c  .  found  in  the  ludcnq)  subroutine,  eg  nmax  should  . 
c  .  have  the  same  numerical  value  as  in  ludcmp 

. . 

complex  *16  a(np,np),alud(np,np),b(n),x(n),sdp,r(nmax) 
dimension  indx(n) 
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do  12  i=l,n 
sdp=-b(i) 
do  1 1  j=l,n 

sdp=sdp+a(i,j)*x(j) 

1 1  continue 
r(i)=sdp 

12  continue 

call  lubksb(alud,n,np,indx,r) 
do  13  i=l,n 
x(i)=x(i)-r(i) 

1 3  continue 
return 
end 


double  precision  function  find_g(nbg,nn,ang,cte,dte,ctm,dtni) 

c  .  This  function  will  calculate  the  asymmetric  parameter  g  which  is  . 
c  .  used  in  determining  radiative  transfer.  It  will  call  the  subrou- 
c  tine  "rotation"  and  use  the  resulting  values  to  calculate  g 
c  Passed  into  this  function  are  the  values  for  nbg,  nn  =  isize, 
c  .  the  incident  angle  ang,  and  the  4  scattering  amplitude  coefficients  . 

. . 

integer  nbg,nn 
double  precision  ang 

complex  *16  cte(0;nn,0:nn),dte(0:nn,0:nn) 
complex  *16  ctm(0:nn,0:nn),dtm(0.nn,0.nn) 
parameter  (size=76) 

. . 

c  .  if  nbg  is  bigger  than  this  size  value, 
c  .  then  change  size  to  equal  at  least  nbg.  . 

. . 

double  precision  tl,t2,fact,fh,asym_g,gg 

complex  *16  a(-size:size,0:size),b(-size:size,0:size),cc 
complex  *16  c(-size:size,0:size),d(-size:size,0:size),ccc 

integer  n,m 

gg  =  0.0 

ccc=dcmpbc(0.0d0, 1 .0) 

ij=size 

c  print*, 'nn,ij-,nn,ij 

call  rotation(a,b,c,d,nbg,ang,cte,dte,ctm,dtm,nn,ij) 

do  10n=0,nbg,l 
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do  20  mp— n,n,  1 
m=mp 

tl=m*dreal(a(iTi,n)*dconjg(b(m,n))  +c(m,n)*dconjg(d(in,n))) 
fti=dfloat(n) 

fact=(  lTi-m+ 1 . )  *(fh+m+ 1 .  )/(2 .  *fo+3 .  )/(2 .  *fh+ 1 . ) 
fact=fh*(fti+2 .  )*dsqrt(fact) 

cc=a(  m,n)*dconjg(a(m,n+ 1  ))+b(m,n)*dconjg(b(m,n+ 1 )) 
cc=cc+c(m,n)*dconjg(c(ni,n+l))+d(m,n)*dconjg(d(m,n+l)) 
t2==fact*dreal(  ccc  *cc ) 
asym_g=4.  *(t  1 +t2) 
gg  =  gg+asym_g 
20  continue 

1 0  continue 

find_g=gg 
return 
end 

subroutine  rotation(a,b,c,d,nbg,ang,cte,dte,ctm,dtni,nn,ij) 

c . 

c  .  This  subroutine  calculates  the  new  rotated  values  of  the  scattering  . 
c  .  coefficients  C's  and  D's,  and  return  these  new  values  to  the  fiiction. 
c  .  find  g  to  find  the  asym  parameter  g. 

c . 

integer  nbg,nn,ij 

complex  *16  a(-ij:ij,0:ij),b(-ij:ij,0:ij) 
complex  *16  c(-ij:ij,0:ij),d(-ij:ij,0:ij) 
double  precision  ang 

complex  *16  cte(0:nn,0:nn),dte(0:nn,0:nn) 
complex  *16  ctm(0:nn,0:nn),dtm(0:nn,0:nn) 

parameter  (siz=200) 

c . 

c  .  the  parameter  siz  here  must  be  at  least  . 
c  .  twice  the  parameter  size  given  in  find  g  . 

c . 

double  precision  factorial 
double  precision  beta(-siz:siz,-siz:siz) 
double  precision  tl,t2,t3,t4,fact,tennl,term2 
complex  *16  suml,sum2,sum3,sum4,cl,c2,c3,c4 
integer  n,mp,m,iijj,i_high,i_low 


c  print  *,'angle=',ang 

ang=-ang 
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do  10  n=0,nbg 
c  print*,'rotation  n=  ',n 
do  15  m=-n,n 
a(m,n)=cmplx(0.0,0.0) 
b(m,n)=cinplx(0.0,0.0) 
c(m,n)=cmpk(0.0,0.0) 
d(m,n)=cnpbc(0.0,0.0) 

1 5  continue 

. . . . 

c  .  For  each  value  of  n,  we  will  find  aU  possible  values 
c  .  of  the  BETA(m,nip),  and  find  the  resulting  new  scat- 
c  .  tering  coefficients  a,b,c,d 

. . 

do  20  mp=-n,n 
do  30  m=-n,n 

1 1  ^actorialfn+mp) 

t2=faetorial(n-mp) 

t3=factorial(n+m) 

t4=factorial(n-m) 

terml=tl+t2-t3-t4 

term2=dexpftennl ) 

beta(mp,m)=dsqrt(term2) 

suml=cnq>bc(0.0,0.0) 

i_low=0 

if((-m-mp).gt.iJow)  iJow=-m-mp 
i_high=n-m 

if((n-mp).lt.i_high)  i_high=n-mp 
if(i_low.gt.i_high)  then 
beta(mp,m)=0.0 
else 

do  40  jj=i_low,i_high,l 
tl=factorial(n+m) 
t2=factorial(n-mp-j[j) 
t3=factorial(m+nip+jj) 
t4=tl-t2-t3 
tennl=dexp(t4) 

tl=factorial(n-m) 

t2=factorial(ij) 

t3=factorial(n-m-j|j) 

t4=tl-t2-t3 

tenn2=dexp(t4) 
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ii=n-m-li 

ii=mod(ii,2) 

if(ii.eq.O)then 

fact=1.0 

else 

fact— 1.0 
endif 

if(terml  .gt.  1  .d200.or.tennl  .It.  1  .d-200)then 

prmt*,'nip,m,  beta=',mp,m,beta(nip,m) 
stop 

endif 

ii=2*jy+mp+m 
t3=(dcosd(ang/2.))**ii 
ii=2*n-2*j[j-mp-m 
t4=(dsind(ang/2.))**ii 
sum  1  =sum  1 +( term  1  *  terni2 )  *fact  *  t3  *t4 
40  continue 

endif 

beta(mp,m)=beta(mp,m)  *  sum  1 

c . 

c  .  That's  one  value  of  m  computed  for  . 

c  .  each  given  value  of  mp 


30  continue 

20  continue 

. . 

c.  There!  We've  just  finished  calculating  all  the 
c.  possible  values  of  BETA(mp.m)  for  each  given  . 
c  .  value  of  n. 


. . 

c  .  Now  to  find  the  new  values  for  the  scattering 
c  .  coefficients,  using  the  BETA  that  we've  just  found 

. . 

do  50  mp=-n,n,l 
suml  =cmplx(0.0,0.0) 
sum2=cmplx(0.0,0.0) 
sum3=cmpbc(0.0,0.0) 
sum4=cn^k(0.0,0.0) 

666  do60m=-n,n,l 

if(m.le.0)  then 
ii=abs(m) 

if(mod(ii,2).eq.O)  then 
cl=cte(ii,n) 
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c2=-dte(ii,n) 

c3=-ctm(ii,n) 

c4=dtm(ii,n) 


else 

cl  =  -cte(ii,n) 
c2=  dte(ii,n) 
c3=  etm(ii,n) 
c4=  -dtm(ii,n) 
endif 
else 

cl=cte(m,n) 

c2=dte(m,n) 

c3=ctm(m,n) 

c4=dtm(m,n) 

endif 

suml  =  suml+cl*beta(m,nip) 
sum2  =  sum2+c2*beta(m,mp) 
sum3  =  suni3+c3*beta(m,mp) 
sum4  =  sum4+c4*beta(m,mp) 

60  continue 

. . 

c  .  These  a.b.c.d  are  now  the  new  values  of  the  scattering  . 
c  .  coefficients  after  the  rotation 
. . 


a(np,n)=suml 
b(mp,n)=sum2 
c(mp,n)=sum3 
d(inp,n)=sum4 
50  continue 
1 0  continue 

return 
end 

^  ********+*****♦****♦******♦♦***♦****•*******♦***** 

double  precision  function  factorial(n) 

integer  n 

integer  loop 

double  precision  sum 

. . 

c.  Since  n  can  be  a  huge  number  (like  190)  we  will  take  . 
c.  the  log  of  the  factorial,  so  the  machine  won't  blow 
c.  up.  eg.  10!  will  be  Iogl0+log9+log8+...logl 
. . 
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sum=0.0 
if  (n.lt.O)  then 

print*, 'ERROR!  value  less  than  zero!' 

print*,'value-,n 

stop 

endif 

if  (n.gt.l)  then 
do  10  loop=l,n 

suni=sum  +  dlog(dfloat(loop)) 

10  continue 

factorial=suin 

else 

factorial=0.0 
endif 
99  return 
end 
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